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COSMIC  RAY  EFFECTS  ON  MICROELECTRONICS 

PART  IV 


1.0  INTRODUCTION 

This  report  is  intended  to  provide  some  of  the  tools  needed  to  estimate 
single  event  upset  (SEU)  rates.  An  SEU  occurs  when  a  single  ionizing 
radiation  event  produces  a  burst  of  hole-electron  pairs  in  a  digital 
microelectronic  circuit  that  is  large  enough  to  cause  the  circuit  to  change 
state.  In  space,  these  unplanned  changes  of  state  usually  result  from  the 
direct  ionization  of  single  heavy  ions  originating  outside  the  spacecraft. 
SEU’s  may  also  be  caused  by  the  products  of  nuclear  reactions  initiated  by 
particles  that  originated  outside  the  spacecraft. 

SEU’s  are  causing  operational  difficulties  of  various  kinds  on  more  than  a 
dozen  spacecraft.  Obviously,  vulnerability  to  SEU’s  must  be  considered  in  the 
design  of  future  spacecraft.  In  order  to  assess  the  vulnerability  of  any 
proposed  design,  engineers  must  have  a  reliable  means  of  estimating  the  SEU 
rates  in  the  radiation  environments  that  can  be  expected  during  the  mission. 
This  requires  experimental  measurements  and  circuit  modeling  to  establish  the 
parameters  that  determine  the  SEU  sensitivity  of  each  device  in  the  circuit. 
It  also  requires  a  reliable  model  of  the  near-earth  particle  radiation 
environment. 

The  first  report  in  this  series  (Adams  et.  al.,  1981a)  provided  a  model  of 
the  near-earth  particle  environment  at  the  orbit  of  the  earth,  but  outside  the 
magnetosphere.  The  second  report  (Adams  et.  al.,  1983)  described  a  method  for 
computing  the  orbit-averaged  geomagnetic  cutoff  transmission  function  that  is 
used  to  modulate  the  environment  outside  the  magnetosphere  to  any  spacecraft 
orbit.  The  third  report  of  the  series  (Tsao  et.  al.,  1984)  described  the 
radiation  environment  high  in  the  earth’s  atmosphere. 

This  report  first  describes  some  improvements  in  the  model  of  the 
near-earth  particle  environment  and  extends  the  method  of  calculating  the 
geomagnetic  cu  off  transmission  function  to  elliptical  orbits.  Next,  the 
method  of  transporting  the  radiation  environment  through  the  walls  of  the 
spacecraft  is  presented.  The  report  goes  on  to  describe  the  transformation 
from  energy  spectra  to  LET  spectra  and  the  computation  of  single  event  upset 
rates  caused  by  the  direct  ionization  of  particles  originating  outside  the 
spacecraft.  Several  examples  of  SEU  rate  calculations  are  given  and  the  SEU 
rates  due  to  direct  ionization  are  compared  with  those  that  result  from  the 
products  of  nuclear  reactions  caused  by  particles  that  originate  outside  the 
spacecraft.  Computer  programs  are  given  that  allow  the  user  to  compute  the 
SEU  rate  of  any  device  inside  any  spacecraft  on  any  orbit. 
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2.0  REVISIONS  AND  EXTENS IONS  OF  PREVIOUS  WORK 

Our  earlier  work  on  cosmic  ray  effects  on  microelectronics  has  been 
revised  and  extended  as  described  in  the  following  sections. 

2.1  Update  of  the  near-earth  particle  environment 

Since  the  publication  of  the  first  report  in  this  series  (Adams  et.  al., 
1981a),  new  data  on  the  near-earth  particle  environment  has  been  published. 
Using  these  results  and  others,  the  model  presented  in  the  first  report  has 
been  improved  and  extended  (see  Appendix  1). 

We  used  the  data  from  the  French-Danish  experiment  on  HEAO-3  (Engelmann 

et.  al.,  1983,  and  Juliusson  et .  al . ,  1983)  to  improve  the  fits  to  the 

galactic  cosmic  ray  data  for  the  elements  up  to  nickel. 

Figures  1  and  2  give  the  new  fits,  to  the  helium  and  iron  differential 

energy  spectra  for  galactic  cosmic  rays.  The  lower  solid  curves  in  each 
figure  are  fits  to  the  pure  galactic  cosmic  ray  spectra  at  solar  maximum. 
This  is  the  mildest  radiation  environment  in  the  11-year  solar  cycle.  As  can 
be  seen  in  figure  2,  no  data  have  yet  been  published  on  the  Fe  spectrum  at 
solar  maximum.  The  upper  solid  curves  are  fits  to  pure  galactic  cosmic  ray 
spectra  at  solar  minimum.  The  galactic  cosmic  ray  flux  is  highest  at  the 
minimum  of  the  solar  activity  cycle.  The  dashed  curves  are  spectra  that 

include  contributions  frcm  co-rotating  interaction  regions,  the  anomalous 
component,  and  small  solar  flares.  These  spectra  are  chosen  to  present  fluxes 
that  are  so  high  at  each  energy  that  flux  levels  exceeding  these  occur  only 
10?  of  the  time.  The  new  fits  should  predict  the  absolute  cosmic  ray  flux 
within  a  factor  of  two.  This  uncertainty  arises  primarily  from  our  inability 
to  predict  future  levels  of  solar  modulation.  The  new  fits  predict  the 
relative  abundances  of  the  elements  below  nickel  to  ±15t,  on  average,  at  the 
energies  where  the  particle  fluxes  are  most  intense. 

The  model  was  extended  to  include  all  the  elements  up  to  uranium.  This 
wa3  done  in  the  following  way:  The  data  from  the  HEAO-3  Heavy  Nuclei 
Experiment  (Binns  et.  al.,  1981;  Binns  et.  al.,  1982,  and  Binns  et.  al.,  1983) 
provided  measurements  of  the  elemental  abundances  for  the  elements  heavier 
than  nickel.  These  measurements  did  not,  however,  resolve  the  individual 
elements.  It  was  necessary  to  determine  the  individual  elemental  abundances, 
in  most  cases,  by  matching  a  cosmic  ray  propagation  calculation  (for  an 
example,  see  Tsao  et.  al.,  1981)  to  the  HEAO  measurements.  The  cosmic  ray 
source  composition  wa3  initally  taken  to  be  the  general  abundance  of  elements 
(Cameron,  1980)  and  then  the  source  composition  and  propagation  conditions 
were  adjusted  to  reproduce  the  measured  abundances  at  earth.  The  ''arriving” 
elemental  abundances  at  earth  from  this  propagation  calculation  were  taken  as 
the  correct  elemental  abundances  for  galactic  cosmic  rays. 

Ihe  description  of  the  contribution  from  the  anomalous  component  (Adams 
et.  al.,  1981)  has  been  extended  to  allow-  investigation  of  the  effect  it  would 
have  on  spacecraft,  were  it  to  be  singly-ionized.  This  has  been  done  by 
independently  modeling  the  anomalous  component  spectra,  based  on  the  data  of 
Webber  et.  al.  (1979).  Both  the  singly-ionized  anomalous  component  and  the 
galactic  cosmic  ray  component  of  each  particle  spectrum  must  be  modulated  to 
the  orb  it- aver aged  geomagnetic  cutoff  (Adams  et.  al.,  1983)  for  the  spacecraft 
according  to  their  different  magnetic  rigidities  at  the  same  energies.  The 
two  modulated  components  are  then  combined  to  form  the  orbit-aversged 
differential  energy  spectrum  incident  on  the  skin  of  the  spacecraft. 


Proton  differential  energy  spectra  for  three  solar  flare  models  were 
presented  in  Adams  et.  al «  (1981a).  The  ordinary  and  worst  case  flare  models 
were  for  large  flares,  i.e.,  those  ^dth  one  week  integral  proton  fluxes  above 
10  MeV  exceeding  2.5x10'  protons/ cm  .  These  models  were  based  on  fits  to  the 
log-normal  distributions  of  fluxes  above  three  energy  thresholds.  We  have 
revised  these  fits  by  extending  the  data  base  for  integral  proton  fluxes  above 
10  MeV  through  April,  1984  using  GOES  satellite  data  (Coffey,  1984).  This 
allowed  the  amplitudes  of  the  model  spectra  to  be  adjusted.  See  eqs.  33  and 
34  of  Appendix  1  for  details. 

The  solar  flare  particle  event  descriptions  in  the  original  model  (Adams 
et.  al.,  1981a)  used  the  proton  spectra,  scaled  by  the  heavy  ion  to  proton 
ratios,  to  obtain  the  spectra  of  heavier  ions.  These  ratios  have  been  revised 
on  the  basis  of  additional  data  on  solar  energetic  particle  composition 
(Mewaldt,  1980;  Zinner,  1980;  Briggs,  1979;  and  Van  Hollebeke,  1975).  They 
have  also  been  extended  to  include  the  spectra  for  all  the  elements  heavier 
than  nickel.  The  mean  composition  was  extended  using  the  data  on  the  general 
abundances  of  the  elements  from  Cameron  (1980).  To  get  the  worst-case 
composition  for  these  same  elements,  we  used  the  charge-dependent  enhancement 
factor  suggested  by  Dietrich  and  Simpson  (1978),  and  assumed  that  it  could  be 
extended  to  uranium.  The  only  data  on  heavy  ion  enrichment  beyond  zinc  comes 
from  tracks  in  meteoritic  and  lunar  olivine  samples  (Goswami  et.  al.,  1980). 
These  data  show  an  even  greater  heavy  ion  enrichment  than  predicted  by  the 
formula  we  used  here.  We  will,  however,  continue  to  use  the  Dietrich  and 
Simpson  formula  until  solar  flare  composition  measurements,  made  with 
detectors  designed  for  this  purpose,  become  available.  In  this  report,  we 
have  dropped  the  descriptions  of  the  mean  flux  from  various  flares  given  in 
Adams  et.  al .  (1981),  because  only  instantaneous  flare  fluxes  are  relevant  to 
the  SEU  problem. 

We  also  introduced  another  flare  spectral  model  (Adams  and  Gelman,  1984). 
This  model  is  intended  to  display  all  the  worst  features  of  all  the  recorded 
flares  in  the  last  thirty  years.  Lingenfelter  and  Hudson  (1980)  investigated 
the  statistical  data  on  solar  flare  particle  fluxes.  They  concluded  that 

there  is  a  maximum  size  for  solar  particle  events  and  that  events  larger  than 
those  thus  far  observed  are  very  rare.  We  therefore  will  assume  that  our 
composite  worst-case  flare  model  is  the  most  extreme  event  that  can  occur.  As 
such,  it  presents  the  kind  of  environment  that  spacecraft  should  be  designed 
to  endure  if  they  are  intended  to  function  without  interruption. 

The  updated  model  is  given  in  Appendix  1.  We  have  compared  this  model 
with  the  results  of  Chenette  and  Dietrich  (1984).  We  find  that  its 

predictions  are  in  agreement  with  the  results  reported  by  these  authors.  This 

updated  model  should  be  compared  with  Appendix  1  of  Adams  et.  al.  (1981a)  to 

see  the  detailed  changes. 

In  recent  years,  evidence  has  been  accumulating  that  solar  energetic  heavy 
ions  may  not  be  fully  ionized.  This  is  certainly  the  case  in  the  0.5  _<  E  £ 
2.5  MeV/u  energy  range,  as  has  been  shown  by  A.  Luhn  et  al.  (1984).  At  the 
higher  energies  of  interest  here  solar  energetic  heavy  ions  may  not  be  fully 
ionized  (as  has  been  generally  assumed  up  to  now).  Although  heavy  ions  with 
energy  >10  MeV/u  would  be  fully  ionized  by  passing  through  sufficient  matter, 
the  available  data  only  place  an  upper  limit  their  pathlength  in  matter 
(Mewaldt  and  Stone,  1983). 
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Fischer  et  al .  (1984)  report  evidence  that  solar  energetic  heavy  ions  in 
the  energy  range  5  <  E  <  20  MeV/u  are  not  fully  ionized.  These  authors  report 
upper  limits  on  the  charge  to  mass  ratio  of  heavy  ions  as  low  as  0.1  (  this 
ratio  is  “0.5  for  fully  ionized  heavy  ions).  Ereneman  and  Stone  (1985)  have 
obtained  indirect  evidence  that  solar  energetic  heavy  ions  in  the  energy  range 
3.5  to  50  MeV/u  have  the  same  distribution  of  charge  states  as  measured  for 
0.5  to  2.5  MeV/u  ions  by  Luhn  et  al.  These  authors  have  shown  that  the 
systematic  abundance  variations  in  SEP  heavy  ions  compared  to  solar  coronal 
abundance  can  be  understood  if  the  change  state  distributions  measured  by  Luhn 
et  al.  are  assumed  for  these  higher  energy  heavy  ions. 

The  SEP  model  described  in  this  report  assumes  that  the  SEP  heavy  ions  are 
fully  ionized.  This  assumption  may  be  incorrect  from  the  evidence  discussed 
above.  If  this  is  the  case,  the  SEU  rates  due  to  SEP's  will  be  systematically 
underestimated  for  spacecraft  in  low  earth  orbit,  because  geomagnetic 
shielding  will  not  be  as  effective  as  the  present  model  assumes. 

Under  the  present  circumstances,  the  charge  state  of  SEP '3  is  uncertain, 
so  it's  not  clear  how  the  model  for  SEP's  should  be  altered  to  account  for  the 
SEP  charge  states.  Therefore,  we  recommend  continuing  to  use  the  present 
model.  A  conservative  calculation  can  always  be  made  by  neglecting  the 
protection  afforded  by  the  geomagnetic  cutoff  (assuming  the  geomagnetic  cutoff 
transmission  function  is  1.0  for  all  energies). 

2.2  Extension  of  the  geomagnetic  cutoff  transmission  function 

The  technique  for  calculating  the  geomagnetic  cutoff  transmission  function 
described  in  Adams  et.  al.  (1983)  has  been  extended  to  elliptical  orbits.  See 
Appendix  2  for  the  details. 

3.0  RADIATION  TRANSPORT  THROUGH  THE  WALLS  CF  THE  SPACECRAFT 

The  orbit-averaged  differential  energy  spectrum  of  each  element  at  the' 
skin  of  the  spacecraft  can  be  found  using  the  model  spectra  from  Appendix  1, 
and  modulating  them  to  the  orbit  of  spacecraft,  using  the  geomagnetic  cutoff 
transmission  function  calculated  by  the  method  described  in  Adams  et.  al. 
(1983)  and  Appendix  2.  To  find  the  differential  energy  spectra  at  the 
microelectronic  components  inside  the  spacecraft,  a  transport  calculation  must 
be  performed.  This  is  best  done  by  the  method  of  Adams  (1983).  This  method 
takes  into  account  the  effects  of  energy  loss  and  particle  losses  through 
total  inelastic  collisions.  No  account  is  taken  of  the  way  in  which  the 
projectile  fragments  from  these  collisions  contribute  to  the  differential 
energy  spectra  of  lighter  ions.  This  omission  results  in  a  systematic 
underestimate  of  the  particle  fluxes,  especially  for  the  elements  from  argon 
through  manganese.  This  underestimate  is  so  small,  it  can  be  neglected  for 
shielding  thicknesses  typically  found  in  spacecraft. 

The  differential  ^nergy  spectrum,  f(E),  inside  the  spacecraft  and  behind  a 
thickness,  t  (in  g/cm  of  aluminum  or  equivalent)  of  shielding  is, 

f (E )  =  f»  (E’)CS(E,)/S(E)]exp(-0t),  and  (1) 

a  z  [5X10-26  n(A1/3  +  271/3  -0.4)2]/27.  (2) 
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Where: 


f'(E')  is  the  differential  energy  spectrum  at  the  skin  of  the  spacecraft; 

E'  is  the  energy  at  the  skin  of  the  spacecraft,  i.e.,  E*  =  R”1  CR(E)  +  . 

tl ,  where  R(E)  is  the  residual  range  of  an  ion  having  an  energy  E  and  R*" 
is  the  inverse  function  of  R(E); 

E  is  the  energy  inside  the  spacecraft; 

S(E)  is  the  stopping  power  of  an  ion  having  an  energy  E; 

A  is  the  atomic  mass  of  the  ion;  and 

n  is  Avogadro's  number. 

The  range-energy  and  stopping  power  data  come  from  Adams  et.  al.  (1987).  The 
tables  needed  for  aluminum  are  reproduced  in  Appendices  3  and  4.  Equations 
(1  )  and  (2)  give  estimates  of  the  differential  energy  spectra  inside 
spacecraft  that  are  satisfactory  for  estimating  SEU  rates  provided  the 
shielding  thickness  does  not  exceed  50  g/cm  .  For  greater  shielding 

thicknesses,  the  method  of  Adams  (1983)  seriously  underestimates  the  SEU  rate. 
For  further  discussion  of  this  point,  see  Adams  (1983).  Precise  estimates  can 
be  made  by  an  exact  transport  calculation,  following  the  methods  of  Tsao  et. 
al.  (1984). 

4.0  CALCULATION  OF  THE  LET  SPECTRUM 

The  'next  step  is  to  calculate  the  linear  energy  transfer  (LET)  spectrum 
from  the  differential  energy  spectra  at  the  microelectronic  components  inside 
the  spacecraft.  LET  is  the  rate  at  which  energy  is  deposited  per  unit 
pathlength  of  an  ionizing  particle.  For  the  present  purpose,  it  is  equivalent 
to  the  rate  of  energy  loss  per  unit  pathlength  for  the  same  particle,  i.e., 
dE/dx  or  stopping  power.  Figure  3  shows  the  stopping  power,  S,  versus  energy 
(taken  from  Adams  et.  al.,  1987)  for  several  ions  in  silicon.  Appendix  5 
gives  these  same  data  in  tabular  form. 

The  transformation  from  a  differential  energy  spectrum  to  a  differential 
LET  spectrum  is  simply, 

f(S)  =  f(E  )CdE/dS],  (3) 

By  examining  Figure  3,  it  can  be  seen  that  Eq.  3  has  three  and  sometimes  four 
singularities,  at  the  points  where  dS/dE  becomes  zero  (and  therefore  dE/dS 
becomes  infinite).  When  the  transformation  is  done  on  a  computer,  the  ratio 
of  finite  differentials  is  used  to  approximate  dE/dS  and  if  the  boundaries  of 
the  intervals  in  E  and  S  are  properly  chosen,  the  ratio  can  be  kept  finite. 

To  see  an  example  of  how  this  is  done,  consider  the  three  iron 
differential  energy  spectra  in  Figure  4.  The  lowest  is  pure  cosmic  ray  iron 
at  solar  maximum,  the  middle  curve  is  pure  cosmic  ray  iron  at  solar  minimum, 
and  the  upper  curve  is  the  90?  worst  case  environment  spectrum  for  iron.  Each 
of  these  spectra  is  behind  0,025  inches  of  aluminum  shielding  and  outside  the 
magnetosphere,  so  there  is  no  geomagnetic  cutoff.  Figure  5  shows  the  result 
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of  the  transformation  of  these  spectra  to  differential  LET  spectra.  Three 
spikes  can  be  seen  in  these  spectra,  corresponding  to  the  three  singularities 
in  dE/dS.  These  singularities  make  differential  LET  spectra  awkward  to  handle 
on  a  computer.  We  prefer,  instead,  to  use  integral  LET  spectra.  Figure  6 
shows  these  same  three  spectra  in  integral  form. 

The  final  step  is  to  repeat  the  application  of  Eq.  3  to  the  differential 
energy  spectra  for  all  the  elements  in  cosmic  rays  (i.e.,  protons  to  uranium) 
and  sum  the  resulting  integral  LET  spectra  to  form  one  composite  integral  LET 
spectrum.  This  spectrum  can  then  be  used  to  estimate  SEU  rates  that  result 
from  the  direct  ionization  of  charged  particles. 

5.0  CALCULATION  OF  SEU  RATES  FROM  THE  DIRECT  IONIZATION  OF  CHARGED  PARTICLES 

An  SEU  occurs  when  a  sufficiently  lacge  burst  of  charge  is  collected  on  a 
critical  node  in  one  of  the  digital  microcircuits  on  a  chip.  The  minimum 
charge  required  to  produce  an  SEU  is  called  the  critical  charge.  This  burst 
of  charge  can  come  from  a  segment  of  the  ionized  trail  left  by  the  passage  of 
an  intensely  ionizing  particle.  It  is  assumed  that  each  critical  node  is 
surrounded  by  a  sensitive  volume  (idealized  as  a  rectangular  parallelepiped) 
and  that  the  charge  deposited  in  this  volume  gets  collected  (Pickel  and 
Blandford,  1980).  The  dimensions  of  the  sensitive  volume  are  related  to  those 
of  the  critical  node.  The  sensitive  volume  is  not,  however,  just  the 
dimensions  of  this  feature.  Charge  is  also  collected  from  the  silicon 
surrounding  the  node  by  diffusion.  The  efficiency  of  charge  collection  from 
beyond  the  node  falls  off  with  distance.  Pickel  and  Elandford  (1  980)  discuss 
how  the  critical  charge  and  the  dimensions  of  the  sensitive  volune  are  found. 
In  general,  experimental  measurements  of  the  operational  SEU  cross  section  for 
the  device,  as  well  as  design  data  supplied  by  the  manufacturer,  are  required. 
It  is  often  necessary  to  interpret  these  data  using  detailed  circuit  modeling 
computer  programs  before  the  device  SEU  parameters  can  be  determined  (see 
Price  et.  al.,  1981,  Zoutandyk,  1983,  and  Zoutendyk  et.  al.,  1984,  for 
example). 

This  simple  model,  described  above,  predicts  that  when  the  critical  charge 
has  been  collected,  an  SEU  will  occur.  The  amount  of  charge  collected  depends 
linearly  on  the  LET  of  the  ionizing  particle  and  the  length  of  its  path  within 
the  sensitive  volume.  There  is,  however,  another  effect  that  extends  the  3ize 
of  the  sensitive  volume.  The  intense  trail  of  ionization  left  by  the  charged 
particle  alters  the  electric  field  pattern  in  the  neighborhood  of  the  feature. 
The  field  forms  a  funnel  along  the  particle*  s  track  and  this  enhances  the 
efficiency  with  which  charge  is  collected  (Hsieh  et.  al.,  1981).  This  funnel 
effect  can  be  partly  accounted  for  in  the  simple  model  discussed  above  if  the 
dimensions  of  the  sensitive  volume  are  experimentally  determined.  Several 
attempts  have  been  made  to  model  the  funnel  effect  explicitly  (see  McLean  and 
Oldham,  1982;  Messenger,  1982;  and  Oldham  and  McLean,  1983).  The  model  of 
Oldham  and  McLean  (1983)  appears  to  correctly  describe  the  experimental  data 
qualitatively,  and  is  in  quanitative  agreement  for  lightly  ionizing  particles. 
This  model  predicts  that  the  collected  charge  will  vary  with  LET  to  the  4/3 
power  and  that  SEU  sensitivity  may  depend  on  the  charge  collection  time,  that 
in  turn  depends  on  LET.  These  are  the  practical  differences  between  the  model 
of  Mclean  and  Oldham  (1983)  and  the  Pickel  and  Blandford  model.  Because 
detailed  models  that  explicitly  include  the  funnel  effect  have  not  yet 
"matured”,  the  simpler  model  of  Pickel  and  Blandford  will  be  used  here. 


ruiwr  < 


This  method  for  estimating  soft  upset  rates  due  to  the  direct  ionization 
by  particles  originating  outside  the  spacecraft  is  given  below  in  the  form 
presented  by  Adams  (1983).  The  upset  rate,  Ng,  in  upsets/bit-sec  is: 

Hnax 

Ng  =  22. 5irA  Q  DC  p(L )  ]F(L)  /L2  dL,  (4) 

22,5Qcrit/pmax 


where , 


A  is  the  surface  area  of  the  sensitive  volume  in  m  , 

Q  is  the  minimum  charge  required  to  produce  an  upset,  in 
pfcocoulombs, 

k  2 

L  =  1.05x10  MeV  cm  /g,  the  highest  LET  any  stopping  ion  can 
defiver, 

2 

P  is  the  largest  diameter  of  the  sensitive  volume  in  g/cm  , 
max 

L  is  LET  in-  MeV  cm2/g, 

2 

F(L)  is  the  integral  LET  spectrum  in  Particles/m  ster-sec, 


and 


D[p(L )]  is  the  differential  pathlength  distribution  in  the  sensitive 
volume  of  each  memory  cell  in  cm  /g, 


where, 


p(L)  »  22.5  Q  ../L  is  the  pathlength  over  which  an  ion  of  LET,  L,  will 
produce  a  charge  Q  cr.  .  The  constant  22.5  is  the  conversion  from  picocoulombs 
to  MeV,  assuming  3?o  ev  per  hole-electron  pair. 

Equation  (4)  contains  the  implicit  assumption  that  the  LET  of  each  ion  is 
essentially  constant  over  the  dimensions  of  the  critical  volume.  Of  course, 
this  is  not  true  for  stopping  ions  very  near  the  end  of  their  range.  Eq.  (4) 
will  assume  that  the  maximum  LET  of  the  stopping  ion  applies  over  its  entire 
pathlength  in  the  sensitive  volume.  This  assumption  can  result  in  calculated 
energy  depositions  that  exceed  the  residual  energy  of  the  ion.  The  problem  is 
especially  acute  for  large  sensitive  volume  dimensions,  and  threshold  LET 
values  just  below  the  maximum  LET  of  an  ion  that  is  much  more  abundant  than 
all  heavier  ions.  Fortunately,  this  circumstance  arises  rarely.  Eq.  (4)  will 
be  accurate  if  the  flux  of  stopping  ions  is  small  compared  to  fast  ions  having 
the  same  LET.  Care  should  be  taken  in  the  use  of  this  formula  when  the 
threshold  LET  is  just  below  the  edge  of  a  "cliff"  in  the  integral  LET  spectrum. 

Eq.  (4)  assumes  one  continuously  sensitive  critical  node  per  bit.  In 
general,  there  may  be  several  critical  nodes  per  bit,  each  with  its  own 
sensitive  volume  dimensions  and  critical  charge.  In  addition,  these  nodes  may 
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only  be  sensitive  part  of  the  time,  making  it  necessary  to  calculate  partial 
upset  rates  for  each  node  and  then  combine  the  results,  weighted  by  the 
fractional  "live"  time  of  each  node.  The  reader  is  referred  to  Nichols  et . 
al .  (1  983  )  for  a  collection  of  test  data  and  to  Pickel  and  Elandford  (1  980), 
Price  et .  al .  (1981),  Kolasinski  et .  al .  (1  983  )*  Zoutendyk  (1  983),  and 
Zoutendyk  et.  al .  (1  984  )  for  examples  of  how  device  SEU  parameters  can  be 
determined . 

6.0  CALCULATION  OF  SEU»S  RESULTING  FROM  NUCLEAR  REACTIONS,  CAUSED  BY  PROTONS 

The  intensely  ionizing  particles  that  cause  an  SEU  can  be  the  fragments  of 
a  silicon  nucleus  struck  by  a  particle  originating  outside  the  spacecraft. 
The  most  common  particle  in  space  capable  of  causing  such  a  nuclear  reaction 
is  the  proton.  At  present,  the  most  practical  method  for  calculating 
pro  to  r>- induced  SEU's  is  to  measure  the  operational  SEU  cross  section  at  one 
proton  energy  and  then  use  the  method  of  Bendel  and  Petersen  (1  983)  to  find 
the  SEU  rate  in  any  proton  environment.  Bende}.  and  Petersen  give  the  SEU 
operational  cross  section  in  upsets  per  proton/ an  per  bit  as, 

0  =  1x1012  (24/A)14  [1-exp(-0.18Y0,5)]4  (5) 

where , 


Y  =  [(18/A)0,5]  (E -A),  (6) 

where  E  and  A  are  in  MeV. 

Proton  differential  energy  spectra  can  be  obtained  from  Sawyer  and  Vette 
(1976)  for  the  trapped  radiation,  and  Appendix  1  of  this  report  for  the 
exomagnetospheric  components.  The  exomagneto  spheric  components  can  then  be 
modulated  to  any  orbit  by  the  method  discussed  in  Adams  et .  al .  (1983)  and 
Appendix  2.  The  combined  proton  differential  energy  spectrum  incident  on  the 
skin  of  the  spacecraft  can  then  be  propagated  to  the  electronics  inside  using 
a  simplified  version  of  Eq.  1,  giver,  below: 

f  (E  )  =  f’  (E  *  )  [S  (E  ’  )/S  (E  )]  (7) 

The  final  step  is  to  integrate  the  cross  section  over  the  proton  differential 
energy  spectrum,  i.e., 


R  s  (1X10“4)  4*  f (E  )o(E  )dE  (8) 


0 

2 

where  f(E)  is  in  protons/(m  ster.sec.  MeV)  and  R  is  in  upsets  per  bit-sec. 

7.0  DESCRIPTION  OF  THE  CREME  COMPUTER  PROGRAMS 

The  CREME  programs  are  a  group  of  FORTRAN  routines  that  calculate 
differential  and  integral  energy  and  LET  spectra  of  cosmic  rays  incident  on 
the  electronics  inside  any  spacecraft  in  any  earth  orbit  and  the  single  event 
upset  rates  that  result.  Input  parameters  for  running  these  programs  describe 
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the  interplanetary  and  magneto spheric  weather  conditions,  the  spacecraft’s 
orbit,  the  shielding  surrounding  the  electronics,  and  the  characteristics  of 
the  device  under  consideration.  Input  data  files  contain  tabulations  of 

stopping  powers  and  ranges  of  cosmic  ray  nuclei  in  aluminum  and  silicon,  and 
geomagnetic  cutoffs.  The  program’s  output  files  containing  energy  and  LET 
spectra,  and  single  event  upset  rates.  Descriptions  and  instructions  on  how 
to  run  the  CREME  programs  are  given  below.  The  various  options  given  by  the 
programs  are  described  here  following  the  descriptions  of  the  programs 

themselves.  Listings  of  all  of  the  CREME  programs  are  given  in  Appendix  7  of 

this  report. 

7.1  How  to  run  the  CREME  programs 

First,  compile  all  the  programs  and  subroutines,  and  link  the  programs 
with  all  appropriate  subroutines.  (See  Sections  7.2  and  7.3  to  see  which 

subroutines  must  be  linked  with  each  program).  For  each  run  of  the  CREME 
programs,  follow  the  below  instructions: 

1.  Decide  all  the  options  you  want  in  your  run.  These  deal  with:  your 
spacecraft’s  orbit;  which  cosmic  ray  environment  to  consider;  material 
shielding  surrounding  your  device;  and  specifics  about  the  device.  (See 
Section  7.5  for  details). 

2.  If  you  wish  to  consider  geomagnetic  cutoff  and  trapped  proton  effects 
(see  Section  7.5),  then  run  the  two  auxiliary  programs  STASS  and  GE0MAG2  (see 
Section  7.3). 

3.  To  get  differential  or  integral  cosmic  ray  flux  for  any  element  versus 
particle  energy,  run  SPEC  (see  Section  7.2). 

4.  For  differential  or  integral  cosmic  ray  flux  versus  particle  LET,  run 
LET  (see  Section  7.2). 

5.  To  find  upset  rates  from  nuclear  reactions  caused  by  protons,  run 
BENDEL  (see  Section  7.2). 

6.  To  find  upset  rates  caused  by  direct  ionization,  you  must  run  two 
programs.  First  run  LET  to  get  an  integral  LET  spectrum;  then  run  UPSET  to 
get  the  upset  rate  (see  Section  7.2). 

7.2  The  main  CREME  programs 

The  main  CREME  programs  are  SPEC,  LET,  BENDEL,  and  UPSET.  SPEC  and  LET 
produce  spectra;  UPSET  and  BENDEL  produce  upset  rates, 

SPEC  outputs  the  differential  and/or  integral  energy  spectra  for  a  single 
element.  SPEC  options  are:  name(s)  of  output  file(s);  number  of  data  points 
in  output  file(s);  the  atomic  number  of  the  ion  whose  spectra  you  want;  the 
year;  whether  to  include  geomagnetic  cutoff  and  trapped  protons;  the 
interplanetary  weather  index;  and  the  thickness  of  material  shielding.  SPEC 
must  be  linked  with  the  subroutines  INSIDE,  CUT,  CRF,  DEDXAL,  and  RAL.  If 
geomagnetic  cutoff  and  trapped  protons  are  included,  the  intermediate  data 
files  GTRANS.DAT  and  STASS.DAT  (created  by  the  auxiliary  programs  GE0MAG2  and 
STASS)  must  be  present.  The  output  files  of  SPEC  are  in  two  columns:  the 
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first  column  is  differential  p  flux  in  particles/(m  .ster.sec) /(MeV/u)  or 
integral  flux  in  particles/(ir  .ster.sec);  the  second  column  is  energy  in 
MeV/u . 

LET  outputs  differential  and/or  integral  LET  spectra  from  a  range  of 
elements.  Unlike  SPEC,  which  simply  tabulates  the  results  of  its  subroutines, 
LET  must  reparameterize  the  fluxes.  Using  stopping  power  functions,  LET 
transforms  the  energy  spectra  into  LET  spectra.  The  LET  program  has  the  same 
options  as  SPEC.  In  addition,  the  range  of  atomic  numbers  of  elements  to  be 
included  in  the  LET  spectre  is  required.  LET  must  be  linked  wiun  the 
subroutines  INSIDE,  CUT,  CRF,  DEDXAL,  RAL,  and  DEDXSI.  As  with  SPEC,  the 
intermediate  data  files  GTRANS.DAT  and  STASS.DAT  are  needed  if  cutoff  and 
trapped  protons  are  to  be  included.  The  output  files  of  LET  have  the  same 
format  as  ttjose  of  SPEC,  but  w^th  different  units:  differential  flux  ig 

perticles/(tn  .ster.secJAMeV.cpi  /g)  ojr  integral  flux  in  particles/( in 
.ster.sec);  versus  LET  in  MeV.cm  /g. 

BENDEL  outputs  the  upset  rate  of  a  microelectronic  device  due  to  nuclear 
reactions  caused  by  protons.  The  program  integrates  the  formula  of  Bendel  and 
Petersen  (1583)  over  the  differential  energy  spectrum  of  protons.  BENDEL 
options  are:  number  of  points  used  to  calculate  the  proton  energy  spectrum; 
the  year;  whether  to  include  geomagnetic  cutoff  and  trapped  protons;  the 
interplanetary  weather  index;  the  thickness  of  material  shielding;  and  either 
Bendel* s  2Parameter  A  or  a  measured  upset  cross-section,  in  upsets/(bit 
proton/ cm  ),  at  some  energy,  in  MeV.  BENDEL  uses  the  same  subroutines  and 
data  files  as  SPEC.  It  outputs  upsets/(bit  .sec)  and  upsets/(bit  .day) . 

UPSET  is  a  more  complicated  program  that  calculates  the  upset  rate  from 
the  direct  ionization  of  particles  originating  outside  the  spacecraft,  UPSET 
assumes  that  each  bit  is  stored  in  a  microcircuit  which  has  one  continuously 
sensitive  node.  This  node  is  surrounded  by  a  sensitive  region  which  is 
assumed  to  have  the  shape  of  a  rectangular  parallelepiped.  To  calculate  the 
upset  rate  using  UPSET,  one  must  first  use  LET  to  create  an  integral  LET 
spectrum.  UPSET  then  integrates  this  spectrum  over  the  differential 
distribution  of  pathlengths  through  the  sensitive  region  of  the  one  sensitive 
node  in  the  microcircuit  to  obtain  the  SEU  rate.  UPSET* s  input  parameters 
are:  the  critical  charge  required  to  produce  an  SEU  in  the  device;  the  linear 
dimensions  of  the  parallel i piped  representing  sensitive  region  of  a 
microcircuit;  and  the  name  of  the  file  containing  the  input  integral  LET 
spectrum.  UPSET  must  be  linked  with  the  subroutine  DIFPLD.  UPSET  outputs 
upsets/(bit  .sec)  and  upsets/(bit  .day) , 

7.3  Auxiliary  CREME  programs 

STASS  provides  input  data  to  LET,  SPEC,  and  BENDEL  on  the  trapped  proton 
flux,  averaged  around  the  spacecrafts'  orbit.  It  is  a  simple  program  that 
transforms  unformatted  hand  input  into  formatted  output  for  use  as  input  data 
in  CREME  main  programs.  Data  for  STASS  should  come  fran  the  "Averaged  Differ. 
Flux,  #/cm**2/sec/keV"  column  of  "Composite  Orbit  Spectrum"  tables.  These  are 
typically  included  in  a  National  Space  Science  Data  Center  report  such  as 
Stassinopoulos  (1981)  or  Stassinopoulos  (1982).  STASS  needs  no  subroutines  or 
input  data  files.  It  outputs  a  two-column  file:  the  first  ^column  is  energy 
in  MeV;  the  second  column  is  differential  flux  in  protons/(cni  sec)/keV.  This 
file  must  be  renamed  to  STASS.DAT  before  it  can  be  used  as  input.  (The 
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subroutine  CRF  converts  the  units  in  this  file  to  the  units  used  in  the  main 
programs).  Note  that  the  trapped  proton  contribution  to  SPEC,  LET,  or  BENDEL 
can  be  omitted  by  creating  a  STASS.DAT  file  with  zero  fluxes. 

GEOMAG 2  is  a  long,  complicated  program  that  computes  the  geomagnetic 
cutoff  transmission  function  corresponding  to  a  given  orbit.  The  program 

traces  the  orbit  in  small  time  steps  for  two  days  of  the  orbit,  calculating 
the  position  and  then  the  cutoff  200  times.  These  cutoffs  are  integrated  into 
a  "transmission  function",  see  Adams  et.  al.  (1983)  and  Appendix  2  of  this 
report  for  details  on  this  process.  GE0MAG2  options  are  just  those  parameters 
of  the  orbit  that  are  mentioned  in  Section  7.5.  They  should  describe  the 
orbit  corresponding  to  the  STASS.DAT  file  being  used.  GE0MAG2  does  not  need 
to  be  linked  with  any  other  program  files;  it  uses  the  input  file  CUT0FF.DAT, 
a  table  of  cutoffs  at  20  km  altitude  (see  Appendix  6).  GE0MAG2  outputs  a 

two-column  data  file:  the  first  column  fs  cutoff  in  GeV/ec;  the  second  column 
is  the  transmission  function,  which  varies  from  0  to  1 .  This  file,  called 

GTRANS.DAT,  is  used  as  input,  along  with  a  corresponding  STASS.DAT  file,  for 

the  main  programs— if  the  "geomagnetic  cutoff  and  trapped  protons"  option  is 
chosen . 

7.4  CREME  subprograms 

INSIDE  is  a  subroutine  that  accounts  for  the  effect  of  material  shielding 
on  particle  fluxes  and  energies.  It  considers  energy  loss  effects  and  nuclear 
fragmentation,  but  does  not  keep  track  of  the  fragments.  The  shielding  is 

assumed  to  be  aluminum.  INSIDE  is  called  by  SPEC,  LET,  and  BENDEL;  it  calls 
CRF,  CUT,  RAL,  and  DEDXAL. 

CUT  applies  the  geomagnetic  cutoff  transmission  function  to  individual 

particle  fluxes  and  adds  the  singly- ionized  anomalous  cosmic  rays  (if  that 
weather  condition  has  been  chosen)  .  CUT  is  only  called  when  the  "geomagnetic 
cutoff  and  trapped  protons"  option  has  been  chosen.  In  this  case,  CUT  is 

called  by  INSIDE;  it  calls  CRF  and  uses  the  GTRANS.DAT  data  file. 

CRF  uses  formulas  to  calculate  the  differential  energy  spectrum  of  solar 
flare  particles  and/or  cosmic  rays  of  a  given  element,  energy,  year,  and 

"weather  condition".  In  addition,  CRF  has  an  entry  point  called  PROTON  that 
extracts  or  bit- aver  aged  trapped  proton  fluxes  from  STASS.DAT.  CRF  is  called 
by  INSIDE  or  CUT;  it  calls  no  subroutines  and  uses  the  data  file  STASS.DAT  (if 
the  "geomagnetic  cutoff  and  trapped  protons"  option  has  been  chosen). 

DEDXSI,  DEDXAL,  and  RAL  are  subroutines  that  interpolate  from  range  and 
stopping  power  tables  (see  Appendixes  3,  4,  5).  DEDXSI  is  called  by  LET; 
DEDXAL  and  RAL  are  called  by  INSIDE.  The  subroutines  use  the  input  data  files 
SILICN.DT1,  ALUMNM.DT1,  and  ALUMNM.DT2,  respectively. 

DIFPLD  computes  the  differential  distribution  of  linear  path  lengths  in  a 
rectangular  parallelepiped,  from  the  formula  of  Bendel  (1984).  DIFPLD  is 
called  by  UPSET;  it  calls  no  subroutines. 
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7.5  Options  and  parameters  in  the  CREME  programs 

The  auxiliary  program  GE0MAG2  requires  an  exact  description  of  the 
satellite's  orbit,  which  the  program  traces.  To  describe  the  orbit,  one  must 
input  the  altitude  at  apogee,  the  altitude  at  perigee,  and  the  orbital 
inclination.  To  describe  the  satellite's  position  at  the  start  of  the 
calculation,  one  must  specify  the  initial  longitude  of  the  ascending  node,  the 
initial  latitude  of  the  ascending  node,  and  the  displacement  of  the  perigee 
from  the  ascending  node.  (This  last  parameter  is  not  needed  for  a  circular 
orbit).  In  addition,  GE0MAG2  asks  whether  there  is  a  "magnetic  storm"  and 
whether  to  take  into  account  the  fact  that  the  earth  casts  a  shadow  on  the 
spacecraft  in  cosmic  rays.  See  Adams,  et.  al.  (1983),  for  details. 

The  main  CREME  programs  (those  that  output  energy  spectra,  LET  spectra,  or 
upset  rates)  require  further  information:  First,  the  cosmic  ray  flux 

subroutines  require  the  date,  the  range  of  nuclear  charges  considered,  and  an 
"interplanetary  weather  index".  The  date  (in  years)  tells  where  one  should  be 
in  the  solar  cycle.  The  weather  index  M  picks  which  of  several  sources  and 
solar  conditions  to  use  in  calculating  the  cosmic  ray  flux  near  the  earth.  A 
weather  index  of  M=1  gives  our  best  approximation  to  the  galactic  cosmic  ray 
flux  at  the  given  date.  This  flux  is  included  in  all  other  "weather 
conditions",  except  for  3.  This  value  (M=3)  gives  worst-case  galactic  cosmic 
ray  fluxes  that  allow  uncertainties  in  flux  data  and  solar  activity.  These 
fluxes  are  so  severe  that  they  have  only  a  10%  chance  of  being  exceeded  by 
actual  fluxes  at  any  moment.  When  M=2,  the  anomalous  component,  assumed  fully 
ionized,  is  added  to  galactic  cosmic  rays.  When  M=4,  a  singly-ionized 
anomalous  component  is  assumed.  The  singly- ionized  particles  are  affected 
differently  than  fully-ionized  particles  by  geomagnetic  cutoff.  When  this 
weather  condition  is  chosen,  the  "geomagnetic  cutoff  and  trapped  protons" 
option  must  be  chosen  as  well. 

Weather  indexes  from  5  through  12  add  solar  flare  particles  to  the 
galactic  cosmic  rays: 

M=5:  peak  ordinary  flare  flux  and  mean  composition; 

Ms6:  peak  ordinary  flare  flux  and  worst-case  composition; 

M=7:  peak  10?  worst-case  flare  flux  and  mean  composition; 

M=8:  peak  10?  worst-case  flare  flux  and  worst-case  composition; 

M=9:  peak  Aug.  4,  1972,  flare  flux  and  mean  composition; 

M=10:  peak  Aug.  4,  1972,  flare  flux  and  worst-case  composition; 

M=11;  peak  composite  worst-case  flare  flux  and  mean  composition; 

Ms12:  peak  composite  worst-case  flare  flux  and  worst-case  composition. 

The  four  flare  options— ordinary,  10?  worst-case,  Aug.  1972,  and 
composite— describe  different  flares.  The  "peak"  flux  is  the  maximum  flux 
resulting  from  the  flare.  "Mean"  and  "worst-case"  composition  describe  flares 
with  normal  and  worst-case  (high)  ratios,  respectively,  of  heavy  ion  fluxes  to 
proton  fluxes.  See  Appendix  1  for  details  on  all  of  these  cases. 

The  thickness  of  shielding  surrounding  the  device  must  be  entered  in  the 
runs  of  the  main  programs.  Shielding  thickness  is  entered  as  the  thickness 
equivalent  in  inches  of  aluminum,  which  the  programs  convert  into  g/cm 
Shielding  of  other  materials  can  be  approximated  by  the  equivalent  thickness 
of  aluninum;  i.e.,  corresponding  to  the  same  thickness  in  mass  per  unit  area. 
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Another  option  of  the  main  programs  is  whether  to  account  for  the 
geomagnetic  cutoff  and  trapped  protons.  This  decision  should  be  made  based  on 
the  spacecraft’s  orbits  calculations  for  low  orbits  should  include  these 
effects,  while  high  orbits  (near  geosynchronous  and  beyond)  need  not.  The 

CREME  programs  do  not  compute  trapped  proton  fluxes}  these  roust  be  entered 

using  the  auxiliary  program  STASS. 

Finally,  the  programs  that  calculate  upset  rates  need  device  parameters. 
BENDEL,  which  computes  the  upset  rate  due  to  proton-induced  fragmentation, 
needs  the  parameter  A.  If  A  is  unknown,  a  measured  value  of  upset 

cross-seotion ,  in  upsets/(bit  proton/ cm; ,  and  the  proton  energy,  in  MeV, 
corresponding  to  this  measurement  may  be  given  (see  Bendel  and  Petersen, 
1983).  UPSET,  which  finds  the  total  upset  rate  due  to  direct  ionization, 
requires  the  linear  dimensions  of  the  rectangular  parallel ipiped  representing 
the  sensitive  region  associated  with  each  bit  and  the  oritical  charge  required 
to  produce  an  upset. 

8.0  DISCUSSION  OF  RESULTS  OBTAINED  WITH  THE  CREME  PROGRAMS 

The  3et  of  programs  described  in  section  7.0  has  been  used  to  investigate 

the  effects  of  interplanetary  weather  conditions  and  geomagnetic  and  material 

shielding  on  SEU  rates.  As  Adams  and  Gelman  (1984)  discussed,  the  critical 
charge  of  a  device  corresponds  approximately  to  a  threshold  value  of  LET 
required  to  produce  an  SEU.  It  follows  then  that  the  particle  flux  above  this 
LET  threshold  is  approximately  proportional  to  the  SEU  rate.  With  this  in 
mind,  integral  LET  spectra  can  be  compared  directly  to  see  their  relative 
effect  on  the  SEU  rate  of  any  proposed  device. 

I  *  '  '  *  '  - 

Adams  (1982)  discusses  how  the  integral  LET  spectrum  changes  under  various 
conditions.  Figure  1  of  his  paper  shows  the  enormous  dynamic  range  in  the  LET 
spectra  that  may  occur  outside  the  earth's  magnetosphere.  The  particle  flux 
behind  0.025  inches  of  aluminum  shielding  and  above  any  LET  threshold  can 
change  up  to  5  orders  of  magnitude  depending  on  the  interplanetary  weather. 
Fortunately,  most  of  this  dynamic  range  is  due  to  infrequent  solar  flare 
particle  events.  The  largest  of  these  occur  only  once  a  decade,  and  are  at 
their  worst  for  only  a  few  hours.  These  facts  lead  to  the  first  important 
conclusion:  "It’s  much  easier  to  build  a  spacecraft  to  work  99%  of  the  time 
than  it  is  to  build  one  that  will  work  100%  of  the  time". 

Figure  2  of  Adams  (1982)  demonstrates  the  effect  of  material  shielding  on 
the  solar  particle  event  resulting  from  the  August  4,  1972  flare  for  the 
elemental  spectra  of  hydrogen  through  nickel.  Figure  2  of  Adams  (1983)  shows 
the  effects  of  material  shielding  on  a  Aug.  4-like  flare  particle  event  with  a 
90%  worst  case  enrichment  in  the  heavy  element  spectra.  These  results  include 
the  contributions  of  the  elemental  spectra  of  all  the  elements  found  in 
nature.  Adams  and  Gelman  (1984)  present  a  composite  worst-case  flare  particle 
event  model.  In  Figure  4  of  their  of  their  paper,  they  show  the  effect  of 
material  shielding  on  this  flare.  In  each  of  the  above  examples,  it  can  be 
seen  that  material  shielding  is  quite  effective  in  reducing  the  SEU  rate  in 
these  severe  environments.  The  second  conclusion:  "Material  shielding  is 
useful  in  moderating  solar  flare  particle  environments  and  reducing  the 
dynamic  range  of  SEU  rates,  where  no  geomagnetic  shielding  is  involved". 


It  should  be  noted,  however,  that  even  large  mounts  of  shielding  cannot 
reduce  the  SEU  rate  to  the  levels  normally  occurring  in  the  interplanetary 
medium.  Comparing  Figures  2  and  3  of  Adams  (1983),  we  can  see  that  4  cm  (1,6 
inches)  of  aluminum  shielding  will  still  leave  the  SEU  rate  from  the  Aug.  4  - 
like  flare  1000  times  worse  than  the  cosmic  ray  background  environment.  If  we 
compare  Figure  3  of  Adams  (1  983  )  with  Figure  4  of  Adams  and  Gelman  (1  984  ),  it 
is  clear  that  even  2  inches  of  aluminum  shielding  will  leave  the  composite 
worst-case  flare  particle  event  10,000  times  more  severe  than  the  cosmic  ray 
background.  Conclusion  1  is  still  valid  for  any  practical  shielding 

thickness. 

Adams  (1  983)  also  investigated  the  effect  of  material  shielding  on  the 
cosmic  ray  background.  Figure  3  of  Adams  (1983)  demonstrates  that  thick 
shielding  is  required  to  obtain  even  modest  reductions  in  the  cosmic  ray 
induced  SEU  rate.  The  third  conclusion:  "Material  shielding  is  ineffective 
against  galactic  cosmic  rays". 

Besides  material  shielding,  the  earth's  magnetic  field  provides  protection 
for  spacecraft  that  are  at  low  altitude  and  near  the  equator.  Adams  (1983) 
shows  how  the  orb it- averaged  LET  spectrum  is  affected  by  orbital  altitude  and 
inclination.  While  SEU's  resulting  from  direct  ionization  will  be  less  for 
lightly  shielded  devices  at  lower  inclinations  over  the  normal  range  of  values 
of  the  critical  charge,  the  degree  of  reduction  will  increase,  generally,  with 
the  critical  charge  (see  Fig.  4  of  Adams,  1983).  Secondly,  at  low 
inclinations,  and  for  lightly  shielded  devices,  there  i3  a  range  of  critical 
charge,  above  2x10  MeV  cm  /g,  for  which  devices  appear  to  be  immune  to  SEU's. 
It  is  unlikely  that  such  devices  would  really  be  free  from  upsets  because  in 
practical  spacecraft  electronic  components  will  not  be  lightly  shielded  from 
all  sides.  Figure  8  of  Adams  0  98:0  demonstrates  that  heavy  shielding  will 

restore  the  LET  spectrum  in  the  2x10  to  1x1(r  MeV  cm  /g  range.  Figure  11  of 
the  same  paper  demonstrates  that  the  SEU  rate  for  a  SBP-9900  in  fact  goes 
through  a  minimum  and  then  increases  at  greater  shielding  depths.  In  such 
cases,  additional  shielding  can  actually  cause  the  SEU  rate  to  increase!  The 
fourth  conclusion:  "One  must  always  consider  the  combined  effects  of  magnetic 
and  material  shielding  on  the  SEU  rate". 

The  effectiveness  of  geomagnetic  shielding  depends  on  the  charge  to  mass 
ratio  of  the  ions.  If  these  heavy  ions  are  completely  stripped  of  their 

electrons,  this  ratio  will  be  approximately  0.5  and  magnetic  shielding  will  be 
as  effective  as  discussed  above.  It  is  certainly  true  that  galactic  cosmic 
rays  are  fully  stripped  of  their  electrons.  The  chemical  and  isotopic 
composition  of  arriving  galactic  cosmic  rays  can  ,only  be  explained  if  these 
particles  have  passed  through  approximately  7  g/cm“  of  interstellar  gas  before 
reaching  us  (Adams  et.  al.,  1981a).  That  is  more  than  enough  matter  to  bring 
these  ions  into  their  equilibrium  charge  state,  and  their  energy  is  so  high 
outside  the  heliosphere  that  their  equilibrium  charge  state  is  fully  ionized. 

There  is  a  component  of  the  particle  radiation  outside  the  magnetosphere 
that  may  be  singly  ionized.  It  is  called  the  anomalous  component.  It  is 

discussed  in  detail  by  Adams  et .  al.  (1981a).  If  this  component  is  singly 

ionized  it  will  penetrate  the  earth's  magnetic  field  with  greater  ease  and 
reach  a  spacecraft  as  lower  energy,  higher  LET  radiation.  The  existing  data 
shows  some  indirect  evidence  that  the  anomalous  component  is  singly  ionized. 
Its  charge  state  will  not  be  known  for  certain  until  an  experiment  currently 
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in  space,  onboard  the  NASA's  LDEF  spacecraft  (Adams  et ,  si.,  1981b)  is 
analyzed.  Figure  11  of  Adams  (1  983)  shows  that  some  devices  like  the  SBP-9900 
are  quite  sensitive  to  the  charge  state  of  the  anomalous  component  when  they 
are  in  a  spacecraft  that's  in  a  28.5°  inclination  orbit.  The  SEU  rate  for  the 
SBP-9900  is  a  couple  of  orders  of  magnitude  higher  than  normal  for  light 
shielding,  if  the  anomalous  component  is  in  fact  singly  ionized.  In  general, 
the  charge  state  of  the  anomalous  component  will  influence  the  SEU  rate  for 
low  altitude  orbits  below  60°  inclination.  The  fifth  conclusion:  ''Until  the 
charge-state  of  the  anomalous  component  has  been  determined,  it's  prudent  to 
assume  that  it's  singly  ionized." 

Unlike  cosmic  rays,  solar  flare  particles  have  not  been  through  much 
matter  before  reaching  earth.  Mewalt  and  Stone  (1983)  aeport  that  solar  flare 
particles  at  50  MeV/u  have  traversed  less  than  30  mg/cnr  of  matter.  In  light 
of  this  result  one  must  consider  the  possibility  that  solar  flare  particles 
are  not  fully  ionized.  Gloeekler  et.  al.  (1981)  measured  the  charge  state  of 
various  solar  flare  ions  from  He  to  Fe  in  the  energy  range  from  0.3  to  2.4 
MeV/u.  The  ionic  charge  states  of  Fe  were  found  to  be  distributed  from  +5  to 
+20,  with  a  broad  maximum  around  +13.  At  higher  energies,  Fischer  et.  al . 
(1984)  recently  reported,  many  ions  in  the  5  to  20  MeV/u  range  from  the  solar 
flare  of  Feb.  13,  1978,  were  not  fully  ionized.  Using  the  geomagnetic  cutoff, 
these  authors  were  able  to  determine  the  maximum  ionic  charge  that  the 
observed  ions  could  have  possessed.  They  report  upper  limits  on  the  ionic 
charge  as  low  as  +3  for  silicon  and  +2  for  oxygen.  Based  on  these  results,  we 
cannot  be  certain  of  the  ionic  charge  of  solar  flare  ions  at  any  energy. 
Further  experiments  will  be  needed  to  establish  the  charge  state  of  solar 
flare  particles  as  a  function  of  atomic  number  and  energy.  The  sixth 
conclusion:  "The  effectiveness  of  geomagnetic  shielding  for  heavy  ions  fnom 
solar  flare  particle  events  is  uncertain.  A  conservative  approach  is  to 
ignore  geomagnetic  shielding  for  solar  flare  heavy  ions." 

Spacecraft  in  earth  orbit  are  exposed  to  the  trapped  radiation  which 
contains  protons  and  possibly  heavier  ions.  The  trapped  proton  population  has 
been  modeled  by  Sawyer  and  Vette  (1  976).  There  may  also  be  significant 
numbers  of  heavier  ions  trapped  in  the  inner  belt.  This  subject  has  been 
explored  by  Adams  and  Partridge  (1982  ).  These  authors  found  that,  based  on 
the  scanty  data  available,  the  contribution  to  the  SEU  rate  from  upsets  caused 
by  trapped  heavy  ions  is  likely  to  be  important  at  least  under  some  conditions 
and  that  they  could  not  rule  out  the  possibility  that  trapped  heavy  ions  may 
cause  SEU  rates  that  are  orders  of  magnitude  larger  than  those  predicted  from 
the  trapped  proton  population.  Our  present  knowledge  of  trapped  heavy  ions  at 
the  energies  needed  to  produce  SEU's  is  so  poor  that  it's  not  possible  to  set 
a  credible  upper  limit  on  them.  Recently,  Adams  et  al .  have  reported 
detecting  trapped  He  ions  but  no  heavier  trapped  ions  at  space  shuttle 
altitudes.  At  the  low  altitudes  investigated  by  these  authors,  the  trapping 
lifetime  for  heavy  ions  is  quite  short  (Blake  and  Fresen,  1977).  As  a  result, 
the  failure  of  Adams  et  al.  to  detect  them  does  not  imply  their  absence  at 
somewhat  higher  altitudes.  Seventh  Conclusion:  "Estimates  of  the  SEU  rate 
for  orbits  below  7000  and  200  m.mi.  are  likely  to  only  be  lower  limits". 

Within  the  magnetosphere,  SEU's  may  result  from  the  direct  ionization  of 
heavy  ions  coming  in  from  outside  the  magnetosphere,  or  from  heavy  ions 
trapped  in  the  earth's  magnetic  field.  SEU's  may  also  result  from  the 
intensely  ionizing  fragments  from  nuclear  reactions  caused  by  protons.  This 
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indirect  mechanism  may  be  important  in  the  inner  Van  Allen  belt  due  to  the 
intense  proton  fluxes  trapped  there.  Unfortunately,  too  little  is  known  about 
trapped  heavy  ions  to  assess  their  relative  contribution  to  SEU  rates  (perhaps 
this  question  will  be  resolved  by  the  Chemical  Release  and  Radiation  Effects 
Satellite  mission).  The  relative  importance  of  the  other  radiation  components 
will  be  investigated  below. 

The  SEU  rates  for  three  devices  will  be  examined  as  a  function  of  orbital 
altitude  for  circular  orbits  at  60°  inclination.  Calculations  were  done  at 
orbital  altitudes  of  200,  400,  600,  800,  1200,  1667,  2593,  3889,  5186,  6389, 
and  10371  km.  The  trapped  proton  environments  were  obtained  from 
StassinopoulosO 981 )  and  Stassinopoulos  and  Bsrth  (1982).  The  parameters  for 
these  devices  are  listed  in  Table  1.  These  parameters  are  intended  only  for 
use  in  the  example  calculations  presented  here  and  may  not  be  accurate  enough 


to  obtain 

exact  SEU  rates. 

* 

Table  1: 

The  Device  Parameters 
8,  9,  11,  12,  and  13. 

Used  to  Calculate 

the  SEU  Rates  in 

Figs.  7, 

DEVICE 

SENSITIVE  VOLUME 

CRITICAL  CHARGE 

"A”  PARAMETER 

REFERENCES 

(in  micrometers) 

(in  picocoulombs) 

HM6508RH 

20x15x2.2 

0.  82 

- 

1 

SBP9989 

10x10x1.8 

0.36  (30%) 

25.2 

2,  3,4 

10x10x1.8 

0.  10(70%x15%) 

5,1 

1  xl  xl.  8 

0.02  (70%x15% ) 

AMD2901B 

80x80x3 

0.25 

23 

5, 1,6,7, 8 

NOTE:  In 

the  SBP9989  entry,  30%  and  70%  refer 

to  the  fraction 

of  critical 

collectors  of  each  type.  The  15/6  refers  to  the  fraction  of  time  the  collector 
is  sensitive.  The  last  two  entries  for  the  SBP9989  reflect  the  reduced  cross 
sections  with  which  these  low  values  of  critical  charge  are  observed. 

REFERENCES: 


1.  Petersen  et.  al.  (1  983). 

2.  Price  et.  al .  (1981). 

3.  Pickel  (1981). 

4.  Pickel  (1983). 

5.  Nichols  et.  al.  (1983). 

6.  Zoutendyk  (1983). 

7.  Zoutendyk  et.  al . ( 1  984  > 

8.  Zoutendyk  (1984). 


Figure  7  compares  the  SEU  rates  for  a  AMD-2901B  from  trapped  protons  and 
heavy  ions  from  outside  the  magnetosphere  as  a  function  of  orbital  altitude 
for  a  60  degree  inclination  orbit.  The  environment  in  the  interplanetary 
medium  is  assumed  to  be  the  90?  worst  case  environment,  and  the  electronic 
components  are  assumed  to  be  shielded  uniformly  by  0.1  inches  of  aluminum. 

Eecause  of  the  large  critical  volume,  this  device  presents  a  large  target  for 


cosmic  rays.  At  the  same  time,  the  relatively  high  critical  charge  limits  the 
sensitivity  of  the  device  to  proton-induced  upsets.  The  result  is  that 
pro  ton- induced  upsets  are  dominant  only  in  the  heart  of  the  inner  belt, 
between  1700  end  3 600  km. 

In  this  calculation,  it  was  found  that  Eq .  (4)  incorrectly  took  into 
account  contributions  fran  trapped  protons.  This  problem  arose,  as  discussed 
in  Section  5,  because  of  the  large  device  dimensions  and  the  value  of  the 
critical  charge.  The  corner  to  corner  distance  across  the  sensitive  volume  is 
113  micrometers.  Multiplying  this  length  by  the  maximum  LET  of  a  stopping 
proton  (see  Appendix  5),  we  get  an  energy  deposition  of  13.8  MeV  or  .62  pc, 
which  is  greater  than  the  critical  charge.  However,  the  energy  of  a  proton 
having  a  residual  range  of  113  micrometers  is  only  3.35  MeV  or  0.15  pC, 
therefore  protons  cannot  actually  cause  SEU's  in  a  AMD2901B  by  their  own 
direct  ionization.  The  problem  was  solved  by  omitting  the  proton  contribution 
from  the  integral  LET  spectra  used  to  "calculate  SEU  rates  with  the  UPSET 
routine. 

Figure  8  gives  the  sare  results,  as  Figure  7  for  the  SBP-9989.  Here  both 
sources  of  SEU's  are  lower,  but  this  time  the  pro  ton- induced  SEU's  are  more 
important.  This  is  because  the  relatively  small  device  dimensions  limit  the 
comsic  ray  contribution,  while  the  low  critical  charge  means  that 
proton- induced  SEU's  depend  more  on  the  device  dimensions.  Figure  9  compares 
the  total  SEU  rates  for  these  devices  with  the  rate  for  the  HM6508RH  that  is 
not  vulnerable  to  proton- induced  SEU's.  The  HM6508RH  is  the  least  sensitive 
of  the  three  devices  considered  here.  Since  the  HM6508RH  is  not  sensitive  to 
pro ton- induced  SEU's,  its  SEU  rate  declines  smoothly  with  orbital  altitude, 
while  the  other  devices  show  peak  SEU  rates  at  about  2600  km.  It  should  be 
emphasized  here  that  Figure  9  should  be  regarded  ss  a  lower  limit  on  the  SEU 
rate,  especially  below  7000  km  since  the  contribution  from  trapped  heavy  ions 
was  not  included. 

So  far  this  discussion  has  delt  with  orbit-averaged  spectra.  It  is 
important  to  note  that  the  geomagnetic  cutoff  is  far  from  constant  around  an 
orbit.  The  geomagnetic  cutoff  will  be  lowest  at  the  northern-  and 
southern-most  extremes  of  the  orbit,  and  highest  at  the  points  where  the  orbit 
crosses  the  geomagnetic  equator.  This  means  that  the  SEU  rate  due  to 
particles  coming  from  outside  the  magnetosphere  will  vary  around  the  orbit, 
more  dramatically  so  at  high  orbital  inclinations.  SEU's  from  particles 
trapped  in  the  magnetosphere  will  also  vary  around  the  orbit.  Stassinopoulos 
(1981)  and  Stassinopoulos  and  Barth  (1982  )  computed  the  fraction  of  mission 
time  thatgsatellites  spend  in  regions  that  are  free  from  trapped  radiation  (<1 
proton/cm  sec).  Figure  10  shows  a  plot  of  one  minus  this  fraction  (i.e.  the 
fraction  of  time  spent  in  the  belts)  ~and  the  fraction  spent  in  the  high 
intensity  regions  of  the  belts  OlxlO13  protons/ an  sec)  as  a  function  of 
orbital  altitude  for  60  degree  inclination  orbits.  From  this  plot,  it's  clear 
that  pro  ton- induced  SEU's  occur  in  only  a  portion  of  the  orbits,  and  a  small 
portion  at  low  altitude  orbits. 

To  explore  this  further,  we  compute  the  instantaneous  SEU  rate  around  one 
orbit  of  a  satellite  in  a  60°  inclination  orbit  at  1667  km  altitude  for  each 
of  the  devices  in  Table  1.  As  above,  the  90%  worst  case  environment  is 
assumed  to  exist  beyond  the  magnetosphere  and  the  components  are  assumed  to  be 
behind  0.1  inches  of  aluminum.  The  instantaneous  trapped  proton  flux  was 
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;en  from  Stassinopoulos  0  981),  Figure  11  compares  the  SEU  rates  due  to 
mic  rays  and  trapped  protons  for  an  AMD-2901B  as  a  function  of  mission 
uime.  We  see  that  trapped  protons  appear  to  dominate  only  in  the  heart  of  the 
radiation  belt  passes.  Cosmic  rays  are  at  their  worst  near  the  magnetic  poles 
and  their  mildest  at  the  equator,  as  expected.  Figure  12  shows  the  same 
comparison  for  the  SBP9989.  Here  the  peak  SEU  rates  in  the  radiation  belt 
passes  are  more  pronounced,  though  lower  than  in  Figure  11.  Figure  13 
compares  the  total  SEU  rates  versus  mission  time  for  the  AMD2901B,  the 
SBP-9989  and  the  HM6508RH.  Because  the  HM6508RH  is  not  vulnerable  to  trapped 
protons,  its  SEU  rate  is  at  a  minimum  when  the  SEU  rates  in  the  other  devices 
are  at  their  peaks.  The  SEU  rates  of  the  AMD2901B  and  the  SBP9989  are  much 
closer  in  the  radiation  belts  than  near  the  magnetic  poles  because  their 
sensitivity  to  cosmic  rays  differs  more  than  their  sensitivity  to 
proton-induced  upsets. 

For  all  these  devices,  however,  the  SEU  rates  are  so  low  that  much  less 
than  one  upset  is  expected  per  orbit  in  each  device.  At  the  SEU  rates 
calculated  for  these  devices,  the  SEU  rate  should  be  no  more  than  one  per 
10,000  bits  per  orbit.  Therefore  variations  in  SEU  rate  around  an  orbit 
cannot  be  easily  observed  in  real  time,  unless  very  sensitive  memory  devices 
are  used  to  construct  very  large  memories.  Conclusion  8:  "While  the 
instantaneous  SEU  rate  can  vary  dramatically  around  an  orbit,  these  variations 
occur  on  time  scales  that  are  usually  small  compared  to  the  mean  time  between 
SEU's.  Therefore  the  orbit-averaged  SEU  rate  is  a  useful  guide  for  designing 
spacecraft  digital  electronics.  Only  when  the  orbit  averaged  SEU  rate  gets  as 
large  as  the  inverse  of  the  orbital  period  should  the  designer  be  concerned 
with  the  instantaneous  rate  around  the  orbit." 
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Figure  1:  The  cosmic  ray  helium  spectrum:  The  solid  curves  are  for  solar 

maximum  (lower)  and  solar  minimum  (upper).  The  dashed  curve  is  the  90% 
worst-case  helium  spectrum.  The  data  points  are  from  Adams  et  al.  (1981a). 
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Figure  5:  The  same  three  cases  as  in  Figure  4,  but  now  the  spectra  have  been 
transformed  into  differential  LET  spectra. 
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ALTITUDE  (km) 

Figure  8:  The  calculated  SEU  rates  for  the  SBP998S.  The  conditions  are  the 
same  as  those  in  Figure  7  and  the  labels  on  the  curves  have  the  same  meanings. 


ALTITUDE  (km) 

Figure  9:  The  total  SEU  rates  (pro ton- induced  +  direct)  for  the  SBP9989  and 
the  AMD2901B  are  plotted  versus  orbital  altitude.  The  conditions  are  the  same 
as  those  for  Figure  7.  This  figure  also  shows  the  total  SEU  rate  versus 
altitude  for  the  HM6508RH  under  the  same  conditions  as  the  other  two  devices. 
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APPENDIX  1:  THE  ANALYTIC  MODEL  FOR  THE  CHARGED  PARTICLE  ENVIRONMENT 

Galactic  cosmic  rays  consist  of  electrons  and  the  nuclei  of  all  the 
elements  in  the  periodic  table;  the  first  28  elements  are  the  most  important 
for  cosmic  ray  effects  on  microelectronics.  These  particles  are  from  outside 
the  solar  system  and  their  flux  at  low  energies  is  anti-correlated  with  solar 
activity  (i.e.  more  cosmic  rays  at  solar  minimum).  The  differential  energy 
spectra  in  particles  per  square  meter  -  sterad|an  -  second  -  million  electron 
volts  per  atomic  mass  unit  (i.e.  particles/(m  .ster  ,sec.(MeV/u) ))  are  given 
in  the  following  paragraphs. 

The  spectra  for  protons  (hydrogen  nuclei),  a-particle  (helium  nuclei),  and 
iron  nuclei  are  given  below  for  energies  above  10  MeV/u: 

•> 

F(E,t)  =  A(E  )  sin[W(t-t0  )]-+B(E),  (1) 

where 

W  =  0.576  radian/ year, 
tQ  =  1  950.  6  A.  D.  date  , 
t  =  current  date  in  years, 

E  =  particle  energy  in  MeV/nucleon, 


IKE)  =  0.5  [f  .  (E)  +  f  v  (E)3,  (2) 

min  max 

A(E )  =  0.5  [f.  (E)  -  f  v  (E)3.  (3) 

min  max 

f  .  and  f  differ  only  by  the  choice  of  constants  in  the  equation 
min  max 

f (E  )  =  10m(E/EQ  )a,  (4) 

where  . 

a  =  aQ  (1  -  exp[-X1  (log1Q  E)°3)  (5) 

and  2 

m  =  C1  exp[-X2  (log-jQ  E)j  -  Cg.  (6) 


The  values  of  the  constants  3q ,  Eg,  b,  X^,  X-,  C1 ,  and  C2  are  given  in  Table  1 
for  each  of  the  elements  hydrogen  (H),  helium  (He) ,  and  iron  (Fe)  for  the 
conditions  of  solar  maximum  and  solar  minimum. 


TABLE  1.  Constants  Used  in  Eqs.  5  and  6  to  Compute  the  Differential  Energy 
Spectra  of  H,  He,  and  Fe  at  Solar  Maximum  and  Solar  Minimum 


Element 

a0 

E0 

b 

X1 

X2 

C1 

C2 

H-min 

-2.20 

1.  1775x10* 

2.685 

0. 117 

0.80 

6.52 

4.00 

H-max 

-2.20 

1.  1775x1  0j 

2.685 

0.079 

0.80 

6.52 

4.00 

He-min 

-2.35 

8. 2700x1 o; 

2.070 

0.241 

0.83 

4.75 

5.  10 

He-mex 

-2.35 

8.2700x1  0* 

2.070 

0.  180 

0.83 

4.75 

5.  10 

Fe-min 

-2.  14 

1.  1750x1 o£ 

2.  640 

0.  140 

0.65 

6.63 

7.69 

Fe-max 

-2.  14 

1.  1750x1 05 

2.  640 

0.  102 

0.65 

6.63 

7.69 
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The  differential  energy  spectra  for  carbon  (C),  oxygen  (0),  fluorine  (F), 
neon  (Ne) ,  sodium  (Na) ,  aluminum  (Al),  and  phosphorus  (P)  are  obtained  by 
multiplying  the  helium  spectrum  (from  Eq .  (1))  by  the  scaling  factors  listed 
in  Table  2. 

TABLE  2.  The  Ratios  of  the  Abundances  of  Various  Nuclei  to  Helium 


Element 

Ratio 

Element 

Ratio 

C 

3.04x10“! 

Na 

1.02x1  0 

0 

2.  84x10"! 

Al 

1.07x1  0 

F 

6.  06x10“Z 

P 

2.  34x1  0 

Ne 

4.  63x1 0"'3 

The  differential  energy  spectra  for  calcium  (Ca),  cobalt  (Co),  and  nickel 
(Ni)  are  obtained  by  mutliplying  the  ihon  spectrum  (from  Eq.  (1))  by  the 
scaling  factors  listed  in  Table  3. 

TABLE  3.  The  Ratios  of  the  Abundances  of  Various  Nuclei  to  Iron 


Element  Ratio 

Ca  2,1x10“! 

Co  3.4x10"z 

Ni  5.0x10** 

The  spectra  of  the  elements  lithium  (Li),  beryllium  (Be),  and  boron  (B) 
are  obtained  from  the  helium  spectrum  FHg  ,  modified  by  the  equation  below. 


0.021  Fu  ,  E  <  3000  MeV/u 

F*  -  He 

-0  44R 

0.729  E  *  3  FHe  .  E  >  3000  MeV/u 


(7) 


F*  is  the  combined  spectrum  of  (Li  +  Be  +  B) .  This  value  is  multiplied  by 
the  ratios  in  Table  4  to  obtain  the  individual  element  spectra. 


TABLE  4.  The  Relative  Fractions  of  Li,  Be,  and  B  in  the  Combined  Total 
Abundance  of  Li  +  Be  +  B 


Element 

Ratio 

Li 

0.330 

Be 

0.176 

B 

0.480 

The  spectrum  of  the  element  nitrogen  (N)  is  obtained  by  modifying  the 
helium  spectrum  F^  as  shown: 


F 


N 


(8.7x10“3  exp[-0.4(log10  E  -  3.152] 

+  7.  6x1  0-3  exp[-0.9(log10  E  -  0.8)2])  FHg. 


(8) 
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The  spectra  of  the  elements  magnesium  (Mg),  silicon  (Si),  and  sulfur  (S ) 
are  obtained  by  modifying  the  helium  spectrum  as  shown: 

*  FHa,  E  >  2200 

F  =  He  .  (9  ) 

(1  +  1.56x10"°  (E-2200))FHe.  E  >  2200 

« 

The  individual  spectra  for  these  elements  are  obtained  by  multiplying  F  by 
the  ratios  of  Table  5. 


TABLE  5. 

The  Ratios  of  Mg,  Si,  and  S  to  an  Adjusted  Helium  Spectrum 

Element 

Ratio 

,  -3 

Mg 

6.02x10  1 

Si 

4.63x10"? 

S 

9.30x10 

The  spectra  for  the  elements  chlorine  (Cl),  argon  (Ar),  potassium  (K), 
scandium  (Sc),  titanium  (Ti) ,  vanadium  (V),  chromium  (Cr),  and  manganese  (Mn) 
are  all  obtained  by  modifying  the  iron  spectrum  Fpe  as  shown  below: 

F*  r  Q  (E) F  ,  where  (10) 

Fe 

Q(E)  *  16  Cl  -  exp(-0.075E0,1*  )]  E"0*33.  (11) 

The  individual  spectra  for  these  elements  are  obtained  by  multiplying  F*  by 
the  ratios  of  Table  6, 

TABLE  6.  The  Fractional  Abundance  of  Each  Element  in  the  Sub-Iron  Group 


Element 

Rat  io 

El ement 

Rat  io 

Cl 

0.070 

Ti 

0. 147 

Ar 

0. 130 

V 

0.070 

K 

0.  090 

Cr 

0. 140 

Sc 

0.042 

Mn 

0.100 

The  differential  energy  spectra  for  elements  from  copper  to  uranium  are 
obtained  by  multiplying  the  iron  spectrum  (from  Eq .  (1))  by  the  scaling 
factors  listed  in  Table  7. 


TABLE  7.  The  Ratio  of  the  Abundances  of  Various  Nuclei  to  Iron 


El ement 

Ratio 

Element 

Ratio 

Cu 

6.8x10  ^ 

Pm 

1.9x10 

Zn 

8.8x10  _ 

Sm 

8.7x10 

Ga 

6.5x10^ 

Eu 

1.5x10 

Ge 

1.4x10 

Gd 

7.0x10 

As 

9.9x10^ 

Tb 

1.7x10 

Se 

5.8x10  J- 

Dy 

7.0x10 

Br 

8.3x10"° 

Ho 

2.6x10' 

Kr 

2.3x1  O’8 

Er 

4.3x1 0“J 
8.9x1  0“? 
4.4x1  0“' 
6.4x1  0’° 

Rb 

1. 1x1  0  “? 

Tm 

Sr 

3.  6x1 O’? 

Yb 

Y 

6.8x10“° 

1 . 7x1  0“? 

Lu 

Zr 

Hf 

4.0x1  0  “' 
3.6x1  O’? 

Nb 

2.6x1  0“? 

Ta 

Mo 

7.  1x10“° 

W 

3.8x1  0“' 
1.3x1 0"' 
5.6x1  0“' 
3.7x1 0“' 
7.2x1 0“I 
1.3x1  o’; 
2.3x1  0“i 

Tc 

1.6x10“° 

Re 

Ru 

5.3x10“° 

0s 

Rh 

1.5x10“° 

Ir 

Pd 

4.5x10"? 

Pt 

Ag 

1.3x1  O’? 

Au 

Cd 

3.6x10“° 

Hg 

In 

1.4x10“° 

T1 

1 . 8x1 0“' 

Sn 

7.5x10“? 

Pb  * 

1.7x1 0’° 
9.0x1 0“° 

Sb 

9.9x1  0“' 

Bi 

Te 

5.7x1  0“? 

Po 

0 

I 

1.5x1 0”? 

At 

0 

Xe 

3.5x10"? 

Rn 

0 

Cs 

5.  8x1 0~l 

Fr 

0 

Ba 

6.0x10“? 

Ra 

0 

La 

5.3x1 0"' 

Ac 

0  „ 

Ce 

1.6x10“? 

Th 

9.0x1  0“° 

Pr 

3.0x1 0“' 

Pa 

0  o 

5.4x1  0’° 

Nd 

1. 1x10“° 

U 

The  recipe  given  shove  is  correct  for  quiet  periods  in  the  interplanetary 
medium  when  only  the  galactic  cosmic  rays  are  present.  These  conditions  are 
often  disturbed,  especially  at  low  energies,  by  small  solar  flares, 
co-rotating  events,  etc.  To  allow  for  typical  disturbed  conditions,  we 
recommend  that  a  worst-case  spectrum  be  used.  With  90  per  cent  confidence, 
the  instantaneous  particle  flux  should  never  be  more  intense  than  described  by 
this  case  at  any  energy. 


To  construct  the  worst-case  spectrum  for  protons,  compute  the 
spectrum  (using  Eq.  4)  and  then  compute  F^-worst  33  shown  bei°WJ 


F 


H -worst 


[l897e"E/9*66  +  1.64  3FU  . 

H— mm 


"H-min" 


(12) 


This  applies  for  E  <  100  MeV/u. 


In  like  manner,  the  solar  minimum  case  helium  and  iron  spectra  (also  from 
Eq.  (4))  are  multiplied  by: 

28.4e“E/13,84  +  1.64  (13) 

for  E  <  100  MeV/u. 


The  worst-case  spectra  of  H,  He,  and  Fe  for  any  element  for  E  >  100  MeV/u 
are  approximated  by  a  multiple  of  the  solar  minimum  spectra: 


F 


worst 


1.64 


Fmin* 


(14) 
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The  resulting  spectra  are  used  as  described  above  to  obtain  the  other 

elemental  spectra,  i.e.  in  the  same  way  as  F„  and  F_  were  used. 

He  re 

In  addition  to  galactic  cosmic  rays,  seme  particles  are  believed  to  be 
accelerated  in  the  interplanetary  medium.  The  most  important  of  these  are 
called  the  anomalous  component.  Its  greatest  contribution  is  to  the  helium 
spectrum.  We  recommend  that,  for  the  period  1  982-1990,  the  cosmic  ray  helium 
spectrum  be  modified  as  follows: 

1.  Determine  the  maximum  values  of  the  helium  spectra  from  Eq.  (4)  using 
the  He-max  and  He-min  constants  from  Table  1 . 


2.  Modify  Eq .  (4)  so  that  each  maximum  value  applies  to  all  energies 
below  the  energy  at  which  the  maximum  occurs;  i.e.,  for  solar  minimum: 


f» 


He-min 


0.33,  E  <  200  MeV/u 

f..  .  (from  Eq.  (4)).  E  >  200  MeV/u 

He-mm 

3.  Make  the  same  kind  of  modification,  f*  ^’e-max  *  *"or  so*ar  maxirnuni: 

0.079,  E  <  300  MeV/u 

He-max  fHe-max  (from  Eq  *  (4))*  E  >  300  MeV/u 

4.  Combine  the  resulting  spectra  as  before  using  Eqs .  (1-3). 


(15) 


(16) 


NOTE:  This  applies  only  to  He;  use  the  regular  He  spectra  of  Eqs.  (1-6) 

to  get  the  spectra  of  the  other  elements. 


The  anomalous  component  contributes  to  the  spectra  of  oxygen  and  nitrogen, 
as  well  as  helium,  at  low  energies.  For  the  years  1982-1  990,  the  following 
contributions  may  be  sdded  to  the  galactic  cosmic  ray  oxygen  and  nitrogen 
spectra . 


For  oxygen,  use: 

f (E )  =  6X10*2  exp'-(ln(E)-1.79  )  2  /0.703  particles/(m2  .ster .sec. MeV/u) .  (17) 


This  spectrum  crosses  the  galactic  spectrum  near  30  MeV/u.  The  two  should  be 
matched  at  that  point  with  Eq .  (17)  replacing  the  galactic  spectrum  at  lower 
energies.  Similarly,  for  nitrogen: 

f(E)  *  1 . 54x1 0-2  exp[g(ln(E  )-1 . 79  )2  /0.70] 

in  particles/ (in  . ster. sec. MeV/u) .  (18) 


Again,  this  crosses  the  galactic  cosmic  ray  spectrum  near  30  MeV/u  and  should 
replace  it  below  this  energy. 

The  spectra  of  the  remaining  elements  are  unaffected  or  affected  at  too 
low  an  energy  to  matter. 


There  is  a  possibility  that  the  anomalous  component  is  singly  ionized.  If 
this  is  so,  it  will  have  an  extraordinary  ability  to  penetrate  the  earth's 
magnetosphere.  In  this  case,  the  differential  energy  speotra  shown  below  are 
assumed  for  the  helium,  carbon,  nitrogen,  oxygen,  neon,  magnesium,  silicon, 
argon,  and  iron  spectra  of  the  anomalous  component.  There  probably  are 
anomalous  components  in  the  spectra  of  the  nuclei  heavier  than  iron,  but  there 
are  no  dat,'  on  them  at  this  time. 


For  singly-ionized  helium, 

.4,  E  <  195 

^  s  4 

1.54x10  E  .  E  >  195 

For  singly- ionized  oarbon, 

«*» 

4.00x103  exp[-(lnE  -  1.79)2  /0.7J,  E  <  10 

F  8  -2 

0.27  E  .  E  >  10 

For  singly- ionized  nitrogen, 

1.54x10"®  exp[-( InE  -  1 . 79 > 2  /0.7],  E  <  20 

F  s  -2 

0.773  E  .  E  >  20 


(19) 


(20) 


(21) 


For  singly- ionized  oxygen, 

6.00x10"®  exp[-(lnE  -  1.79)2  /0.73, 

F  3  -2 

1.32  E  . 

For  singly- ionized  neon, 

S.OOxlO-3  exp(-(lnE  -  1.79)2  /0.7], 

F  3  -2 

0.40  E^  . 

For  singly- ionized  magnesium, 

8.00X10"4  expC-(lnE  -  2.30®  /0.73, 

F  3  -2 

0.16  E  . 

For  singly- ionized  silicon, 

I.OOxlO"3  exp  [-(InE  -  2.20)2  /0.4], 

F  3  -2 

0.10  E  . 

For  singly- ionized  argon, 

5.40X10"4  exp  [-(InE  -  1 . 79 ) 2  /0.7], 


E  <  30 
E  >  30 

E  <  20 
E  >  20 

E  <  20 
E  >  20 

E  <  10 
E  >  10 

E  <  20 


(22) 


(23) 


(24) 


(25) 
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F  s 

(26) 

0.028  E  . 

E  >  20 

For 

singly- ionized  iron, 

6.00X10"4  exp[-(lnE  -  2.48)2  /2.03, 

E  <  30 

F  = 

(27) 

0.35  E^  . 

E  >  30 

Since 

these  anomalous  component  particles 

are  assumed 

to  be 

singly-ionized,  they  will  have  a  higher  magnetic  rigidity  than  galaotic  cosmic 
rays  of  the  same  energy.  The  magnetic  rigidity  of  galaotic  cosmio  rays  is 

R  =  (A/Z)  (E2  +  1862.324  E)1/2  /1 000,  (28) 

in  GeV/ec.  The  rigidity  of  singly  ionized  nuclei  is 

R  s  A  (E2  +  1862.324  E)1  /2  /1 000,  (29) 

To  add  the  singly- ionized  anomalous  component,  it  is  necessary  to  modulate 
both  the  galactic  cosmic  ray  spectra  and  the  anomalous  component  spectra  given 
above  using  the  geomagnetic  cutoff  transmission  function.  The  method  for 
calculating  this  function  is  discussed  in  Adams  et.  al .  (1983)  and  in  Appendix 
2  of  this  report.  The  resulting  modulated  spectra  are  then  added  together. 

Solar  flare  particle  events  are  sporadic  occurrences  lasting  1-5  days. 
When  these  events  occur  they  can  be  the  dominant  oause  of  SEU's.  For 
statistical  treatment,  they  are  broken  down  into  two  classes}  ordinary  (OR) 
and  anomalously  large  (AL).  The  probability  of  having  more  than  a  number  n  of 
events  in  a  time  t  is  given  by: 

P(n,t,N,T)  s  1  .  I  (i4W)!(t/T)i  /[i!N!(1  +  t/T)1+i+N  3,  (30) 

i=0 

where  T  and  t  are  in  years,  and  N  is  the  number  of  flares  that  have  occurred 
in  T  years. 

For  ordinary  events,  Eq.  (30)  becomes: 

Pnp  =  P(n,t,24,7)  for  1981-1983 

and  0R  (31) 

PQR  =  P(n  ,t,  6,8)  for  1984-1988, 

where  there  is  a  probability  P^  of  having  more  than  n  ordinary  events  in  t 
years.  Similarly,  for  anomalously  large  events: 

PAL  =  P(n,t,1,7).  (32) 

If  there  is  an  unacceptable  risk  of  an  AL  event  then  it  will  be  the  worst-case 
flare  for  the  mission. 


The  peak  proton  flux  differential  energy  speotrum  for  ordinary  events  is, 
typically 


FqR  s  2.45x10^  (e"*^2^*^  +  173e"^^  )  protons/m2  .ster .seo .MeV. 

and  no  worse  than 


(33) 


*W0R  *  2.06x10^  (e"^2***^  +  63. )  protons/m2  .ster .seo.MeV,  (34) 

with  a  confidence  of  about  90  per  cent. 

Using  the  August  1972  flare  as  a  model  A1  event,  the  peak  proton  flux 
differential  energy  speotrum  is 


9.3x10*  (dP/dE)  exp(-P/0. 10),  E  <  150  MeV 


'AL  1.76 x105  (dP/dE)  P“9, 


(35) 


E  >  150  MeV 


in  protons/m  .ster .sec  .MeV,  where 

P  *  [(E/1000)2  +  1.86  x  10“3  E]1/2. 


(36) 


To  model  the  worst  flare  that  is  ever  likely  to  occur,  we  use  the 
composite  of  the  two  worst  flares  ever  recorded:  the  August,  1972,  flare  and 
the  February,  1956,  flare.  The  composite  worst-c8se  flare  proton  speotrum  is 
taken  to  be  the  peak  of  the  1972  speotrum,  as  given  by  Eq.  (35),  and  the  1956 
peak  spectrum,  given  below. 

f1956  =  1*1l6x1°8  (E”1*248  )  (0. 248  +  2.5x105  xl.TxEPOW)  EXPOW 

+  4.7x1019  (E"5*3  ) (4.3 x(1 -EXPOW)  -  6.32x10"15  x4.7xEP0WxEXP0W)  ,  (37) 

where  .  - 

EPOW  =  E 

and  c 

EXPOW  =  exp(-2.5  xl  0"°  EPOW). 

The  composition  of  flare  particles  also  varies  greatly  from  flare  to 
flare.  Table  8  gives  the  composition  relative  to  hydrogen  of  the  elements 
through  nickel.  Both  mean  and  (90  percent  confidence  level)  worst-cases  are 
given.  Multiply  the  abundance  ratio  from  Table  8  by  the  appropriate  flare 
proton  spectrum  to  get  the  flare  spectrum  of  any  element  in  the  table. 

TABLE  8.  Mean  and  Worst-Case  Flare  Compositions 


Mean  Case 

Worst  Case 

Mean  Case 

Worst  Case 

H 

He 

I.OxlO"2 

3.3X10"2 

p 

s 

2.3  x10”I 
8.0x10"! 

1.1x10"^ 

5.0x10"! 

Li 

0 

0 

Cl 

1.7x10": 

8.0x10“,: 

Be 

0 

0 

Ar 

3.3x10"! 

1.8x10"! 

B 

0  -4 

0  a 

K 

1.3x10": 

6.0x10“,: 

C 

1.6x10  - 

4.0x10"7 

1.1x10, 

Ca 

3.2x10“° 

2.0x10“3 

N 

3.8x10  jj 
3.2x10“* 

Sc 

0  7 

0  7 

0 

1.0x10“^ 

Ti 

1.0x10  1 

5.0x10  ' 
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Ne 

5.1X10”8 

u  -4 
•  1,9x10  Z 

V 

Cr 

5,7  xIO’Z 

4.0x10 

Na 

3.2x10*2 

1.3x10*2 

2.5x10*2 

Mn 

4.2x10’,! 

2.3x10 

Mg 

6.4x10l 

Fe 

4.1x10’; 

4.0x10 

A1 

3.5x10*2 

1.4x10*2 

1.9x10”* 

Co 

1.0x10*2 

5.5x10 

Si 

5.8x10"° 

Hi 

2.2x10“° 

2.0x10 

The  mean  case  compositions  for  the  elements  from  copper  to  uranium  are 
taken  from  Cameron  (1980).  The  ratios  of  these  abundances  to  hydrogen  are 
given  in  Table  9. 


TABLE  9.  Mean  Flare  Compositions 


Cu 

Zn 

Ga 

Ge 

As 

Se 

Br 

Kr 

Rb 

Sr 

Y 

Zr 

Nb 

Mo 

Tc 

Ru 

Rh 

Pd 

Ag 

Cd 

In 

Sn 

Sb 

Te 

I 

Xe 

Cs 

Ba 

La 

Ce 

Pr 

Nd 


2.  0x1  O' 
6.0x10 
2.0x10 
5.0x1 0 
3.0x10 
3.0x10 
4.0x10 
2.0x10 
3.0x10 
1 . 0x1  0 
2.0x10 
5.0x10 
4.0x10 
2.0x10 
0 

9.0x10 

2.0x10 

6.0x10 

2.0x10 


7.0x10 

9.0x10 

2.0x10 

1.4x10 

3.0x10 

6.0x10 

2.7x10 

2.0x10 

2.0x10 

2.0x10 

5.0x10 

8.0x10 

4.0x10 


-8 

-8 

-9 

-9 

-10 

-9 

-10 
-9 
-10 
-9 
-10 
—1 0 
-11 
-10 

-11 

-11 

-11 

-11 

-11 

-12 

-10 

-11 

-10 

-11 

-10 

-11 

-10 

-11 

-11 

-12 

-11 


Pm 

Sm 

Eu 

Gd 

Tb 

Dy 

Ho 

Er 

Tm 

Yb 

Lu 

Hf 

Ta 

W 

Re 

0s 

Ir 

Pt 

Au 

Hg 

T1 

Pb 

Bi 

Po 

At 

Rn 

Fr 

Ra 

Ac 

Th 

Pa 

U 


■  1.0x10 
4.0x10 
2.0x10 
3.0x10 
2.0x10 
4.0x10 
1.0x10 
2.0x10 
9.0x10 
2.0x10 
8.0x10 
9.0x10 
1.0x10 
2.0x10 
3.0x10 
3.0x10 
6.0x10 
1.0x10 
1.0x10 
9.0x10 
1.0x10 
6.0x10 
0 
0 
0 
0 

..  0 

0 

2.0x10 
0 


-11 
-12 
-11 
-12 
-11 
-12 
-11 
-12 
-12 
-12 
-12 
-13 
—1  1 
-12 
-11 
-11 
-11 
-11 
-11 
-12 
-10 
-12 


-12 


1.2x1 0 


-12 


The  worst-case  compositions  of  the  elements  from  copper  to  uranium  are 
obtained  by  multiplying  the  abundance  ratios  of  Table  9  by: 

(Cy  (0)/Ct  (O  )  )0. 48  exp<Z0,78  /6.89),  (38) 

where  C„  (O)  and  C  (0)  are  the  worst-case  and  mean  abundance  coefficients  for 
oxygen  in  Table  8. 


There  ere  several  good  mathematical  models  for  the  trapped  proton 
environment.  We  reoommend  "AP-8  Trapped  Proton  Environment  for  Solar  Maximum 
and  Solar  Minimum  Donald  M.  Sawyer  and  James  I.  Vette,  Report  No. 

NSSCS/WDC-4-R  and  S  76-06,  Dec.,  1976,  NASA-Goddard ,  Greenbelt,  Md . 

The  environment  of  trapped  a  particles  and  heavier  nuolei  is  very  poorly 
known.  The  reader  is  referred  to  Adams  and  Partridge  (1982)  for  a  discussion 
of  the  data  base  on  trapped  ions  heavier  than  protons. 

The  modulation  of  cosmic  ray  spectra  by  the  earth's  magnetic  field  is 

discussed  in  Adams  et.  al.  (1983)  in  Appendix  2  of  this  report.  Briefly,  the 

geomagnetic  cutoff  is  a  value  of  magnetic  rigidity  below  which  cosmic  rays 

will  not  reaoh  a  specified  point  in  the  magnetosphere  from  a  specified 
direction.  The  magnetic  rigidity  P  in  GeV/ec  may  be  computed  from  a 
particle's  energy  by: 

P  =  (A/Z)[  (E/1000)2  +  1.86x10“3  E]1/2  ,  (39) 


where  A  and  Z  are  the  atomic  mess  and  charge  of  the  nucleus  in  question. 


The  cutoff  at  any  point  for  particles  arriving  from  the  zenith  is  most 
simply  computed  with: 

P  =  15.96/L2*005  GeV/ec,  (40) 

c 

where  L  is  Mcllwain's  L  parameter  (i.e.  the  radial  distance,  in  earth  radii, 
from  the  center  of  the  earth  to  the  point  in  the  geomagnetic  equitorial  plane 
where  it  is  crossed  by  the  magnetic  field  line  that  also  passes  through  the 
point  of  observation) .  .  - 

4  * 

Detailed  calculations  of  the  cutoff  are  available  from  Shea  and*  Smart 
(1975).  Transmission  functions  for  satellite  orbits  may  be  computed  using  the 
techniques  described  in  Adams  et.  al.  (1983)  and  Appendix  2. 


The  transmission  functions  are  useful  in  modulating  the  cosmic  ray 

spectra.  Some  thought  must  be  given  to  their  use  on  solar  flare  spectra 

because:  (1)  the  flare  particle  intensity  changes  on  a  time  scale  comparable 

to  or  shorter  than  an  orbital  period?  (2)  there  is  no  certain  proof  that  solar 

flare  particles  are  fully  ionized;  (3)  the  geomagnetic  cutoff  is  suppressed  to 

some  extent  during  a  flare.  We  recommend  that  the  geomagnetic  cutoff  during  a 

flare,  Pe  ,  be  computed  from  the  "quiet  time"  cutoff  Prt  using: 
r  o 


and 


AP/P  =  0.54  exp(-P  /2.9) 
o  o 


-  A  P, 


where  Pp,  PQ,  and  AP  are  in  GeV/ec. 


(41) 


(42) 


APPENDIX  2:  THE  GEOMAGNETIC  CUTOFF  TRANSMISSION  FUNCTION 

The  geomagnetic  cutoff  transmission  function  is  calculated  using  the 
method  described  in  Adams  et.  al.  (1983).  Here  we  will  describe  how  this 
method  has  been  extended  to  elliptical  orbits. 

Section  4.0  of  Adams  et.  al.  (1983)  describes  the  procedure  for  finding 
the  geographic  position  of  satellites  in  circular  orbits.  For  spacecraft  in 
elliptical  orbits,  the  offset  dipole  coordinates  of  the  position  must  be 
calculated  in  three  steps.  First,  the  position  of  the  satellite  along  its 
orbit  must  be  calculated  versus  time.  Second,  the  geographic  latitude  and 
longitude  of  the  sub-satellite  point  must  be  found.  Third,  the  geographic 
coordinates  must  be  transformed  into  the  offset  and  the  tilted  dipole 
coordinates  required  by  the  Stormer  approximation. 

Numerous  celestial  mechanic  texts  describe  elliptical  orbits  of  satellites 
constrained  by  Newton’s  laws.  For  convenience,  this  description  will  follow 
the  notation  of  Sterne  (I960).  Elliptical  positions  are  described  by  an  angle 
n  called  the  "true  anomaly"  and  the  distance  r  from  the  satellite  to  the 
earth's  center.  These  in  turn  are  formulated  in  terms  of  the  "mean  anomaly" 
M,  the  "eccentric  anomaly"  E,  the  orbit's  period  P,  eccentricity  e,  and 
semi-major  axis  a. 

The  period  is: 

P  =  2ir  (a3  /CM  )1/2, 
e 

where  G  is  the  gravitational  constant  and  M  is  the  earth's  mass.  The  mean 
anomaly  is  given  by:  e 

M  =  (2ir/P )  t, 

where  t  is  the  time  since  perihelion.  The  eccentric  anomaly  is  then  found  by 
iteratively  solving  the  equation: 

M  s  E  -  e  sin  E, 

The  values  r  and  a  are  calculated  from  E: 

r  =  a  (1  -  ecosE), 

tan  (n/2)  *  ((1+e)/(1-e))1/2  tan  (E/2). 

This  last  equation  is  solved  uniquely  with  the  added  condition  that  n/2  and 
E/2  be  in  the  same  quadrant.  For  purposes  of  calculating  geographic 
coordinates,  the  displacement  of  the  perigee  from  the  ascending  node  is  added 
to  this  to  obtain  a  final  value  of  n. 

The  spacecraft's  geographic  latitude  e  and  longitude  $  are  calculated  from 
the  above  r  and  n,  and  the  parameters  of  the  orbit:  the  orbital  inclination 
9q  {  the  initial  longitude,  ,  of  the  ascending  node,  and  the  angular 
velocity,  w,  of  the  earth.  The  latitude  is  simply: 

e  =  ir/2  -  cos”1  (sin©0  sinn). 
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The  longitude  is: 


where 


and 


*  -  -  tf/2  -  wt  +  tan”1  S1  +  tan’1  $2 

sintl  /2  (it/2"-9q  )]cosC1  /2  (ir/2-(j)] 

1  sintl  /2(ir/2+e0)]sinCl/2(tt/2-n)]’ 
cos[1/2(1T/2-e0)3cosCl/2(Tr/2-fi)] 


costl  /2  (,/2+Sq  )] sintl  /2  (ff/2-n)3 

The  final  step — transforming  to  offset-tilted 
using  the  method  of  Smart  and  Shea  (1977). 


dipole  coordinates— is  done 


Appendix  B  of  Adams  et.  al.  (1983)  contained  some  typographical  errors. 
This  appendix  is  presented  in  its  correct  form  as  Appendix  6  of  this  report. 
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APPENDIX  6:  THE  WORLDWIDE  VERTICAL  GEOMAGNETIC  CUTOFF  AT  20KM  ALTITUDE 


LATITUDE 

LONGITUDE 

RIGIDITY 

LATITUDE 

LONGITUDE 

RIGIDITY 

(GeV/eo) 

(GeV/ec) 

80.00 

0.00 

0.02 

70,00 

0.00 

0.26 

80.00 

15.00 

0.04 

70,00 

15.00 

0.34 

80.00 

30.00 

0.06 

70.00 

30.00 

0.41 

80.00 

45.00 

0.09 

70.00 

45.00 

0.47 

80.00 

60.00 

0.09 

70.00 

60,00 

0.49 

80.00 

75.00 

0.10 

70.00 

75.00 

0.51 

80.00 

90.00 

0.10 

70.00 

90.00 

0.52 

80.00 

105.00 

0.11 

70.00 

105.00 

0.55 

80.00 

120,00 

0.11 

.  70.00 

120.00 

0.59 

80.00 

135.00 

0.11 

70.00 

135.00 

0.60 

80.00 

150.00 

0.08 

70.00 

150.00 

0.62 

80.00 

165.00 

0.06 

70.00 

165.00 

0.56 

80.00 

180.00 

0.04 

70.00 

180.00 

0.47 

80.00 

195.00 

0.00 

70.00 

195.00 

0.36 

80.00 

210.00 

0.00 

70.00 

210.00 

0.23 

80.00 

225.00 

0.00 

70.00 

225.00 

0.13 

80.00 

240.00 

0.00 

70.00 

240.00 

0.06 

80.00 

255.00 

0.00 

70.00 

255.00 

0.00 

80.00 

270.00 

0.00 

70.00 

270.00 

0.00 

80.00 

285.00 

0.00 

70.00 

285.00 

0.00 

80.00 

300.00 

0.00 

70.00 

300.00 

0.00 

80.00 

315.00 

0.00 

70.00 

315.00 

0.05 

80.00 

330.00 

0.00 

70.00 

330.00 

0.11 

80.00 

345.00 

0.00 

70.00 

345.00 

0.18 

75.00 

0.00 

0.10 

65.00 

0.00 

0.58 

75.00 

15.00 

0.14 

65.00 

15.00 

0.72 

75.00 

30.00 

0.18 

65.00 

30.00 

0.80 

75.00 

45.00 

0.20 

65.00 

45.00 

0.89 

75.00 

60.00 

0.23 

65.00 

60.00 

0.93 

75.00 

75.00 

0.25 

65.00 

75.00 

0.97 

75.00 

90.00 

0.25 

65.00 

90.00 

1.01 

75.00 

105.00 

0.26 

65.00 

105.00 

1.03 

75.00 

120.00 

0.27 

65.00 

120.00 

1.12 

75.00 

135.00 

0.28 

65.00 

135.00 

1.19 

75.00 

150.00 

0.26 

65. 00 

150.00 

1.20 

75.00 

165.00 

0.24 

65.00 

165.00 

1.13 

75.00 

180.00 

0.20 

65.00 

180.00 

0.95 

75.00 

195.00 

0.14 

65.00 

195.00 

0.74 

75.00 

210.00 

0.09 

65.00 

210.00 

0.53 

75.00 

225.00 

0,03 

65.00 

225.00 

0.32 

75.00 

240.00 

0.00 

65.00 

240.00 

0.17 

75.00 

255.00 

0.00 

65.00 

255.00 

0.09 

75.00 

270.00 

0.00 

65.00 

270.00 

0.05 

75.00 

285.00 

0.00 

65.00 

285.00 

0.04 

75.00 

300.00 

0.00 

65.00 

300.00 

0.08 

75.00 

315.00 

0.00 

65.00 

315.00 

0. 16 

75.00 

330.00 

0.02 

65.00 

330.00 

0.28 

75.00 

345.00 

0.07 

65.00 

345.00 

0.42 

LATITUDE  LONGITUDE  RIGIDITY 

(GeV/ec) 


LATITUDE  LONGITUDE  RIGIDITY 

(GeV/ec) 


60.00 

0.00 

1.14 

50.00 

0.00 

3.21 

60.00 

15.00 

1.34 

50.00 

15.00 

3.54 

60.00 

30.00 

1.46 

50.00 

30.00 

3.81 

60.00 

45.00 

1.57 

50.00 

45.00 

3.97 

60.00 

60.00 

1.61 

50.00 

60.00 

4.14 

60.00 

75.00 

1.67 

50.00 

75.00 

4.27 

60.00 

90.00 

1.73 

50.00 

90.00 

4.36 

60.00 

105.00 

1.82 

50.00 

105.00 

4.37 

60.00 

120.00 

1.95 

50.00 

120.00 

4.68 

60.00 

135.00 

2.05 

50.00 

135.00 

4.93 

60.00 

150.00 

2.05 

*  50.00 

150.00 

4.92 

60.00 

165.00 

1.99 

50.00 

165.00 

4.67 

60.00 

180.00 

1.75 

50.00 

180.00 

4.27 

60.00 

195.00 

1.40 

50.00 

195.00 

3.38 

60.00 

210.00 

1.00 

50.00 

210.00 

2.81 

60.00 

225.00 

0,65 

50.00 

225.00 

2.03 

60.00 

240.00 

0,40 

50.00 

240.00 

1.41 

60.00 

255.00 

0,22 

50.00 

255.00 

0.95 

60.00 

270.00 

0.16 

50.00 

270.00 

0.73 

60.00 

285.00 

0.14 

50.00 

285.00 

0.69 

60.00 

300.00 

0.21 

50.00 

300.00 

0.89 

60.00 

315.00 

0.38 

50.00 

315.00 

1.34 

60.00 

330.00 

0.59 

50.00 

330.00 

1.98 

60.00 

345.00 

0.86 

50.00 

345.00 

2.65 

55.00 

0.00 

1.94 

45.00 

0.00 

4.77 

55.00 

15-00 

2.28 

45.00 

15.00 

5.12 

55.00 

30.00 

2.47 

45.00 

30.00 

5.36 

55.00 

45.00 

2.61 

45.00 

45.00 

5.51 

55.00 

60.00 

2.68 

45.00 

60.00 

5.73 

55.00 

75.00 

2.78 

45.00 

75.00 

5.90 

55.00 

90.00 

2.85 

45.00 

90.00 

6.11 

55.00 

105.00 

2.92 

45.00 

105.00 

6.29 

55.00 

120.00 

3.12 

45.00 

120.00 

6.57 

55.00 

135.00 

3.31 

45.00 

135.00 

6.86 

55.00 

150.00 

3.35 

45.00 

150.00 

6.86 

55.00 

165.00 

3.15 

45.00 

165.00 

6.33 

55.00 

180.00 

2,88 

45.00 

180.00 

5.59 

'55.00 

195.00 

2.22 

45.00 

195.00 

4.85 

55.00 

210.00 

1,75 

45.00 

210.00 

4.08 

55.00 

225.00 

1.23 

45.00 

225.00 

3.16 

55.00 

240.00 

0.78 

45.00 

240.00 

2.37 

55.00 

255.00 

0,50 

45.00 

255.00 

1.74 

55.00 

270.00 

0.36 

45.00 

270.00 

1.32 

55.00 

285.00 

0.36 

45.00 

285.00 

1.22 

55.00 

300.00 

0.46 

45.00 

300.00 

1.49 

55.00 

315.00 

0,75 

45.00 

315.00 

2.21 

55.00 

330.00 

1.13 

45.00 

330.00 

3.16 

55.00 

345.00 

1.59 

45.00 

345.00 

4.20 

LATITUDE 


40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

40.00 

35.00 

35.00 

35.00 

35.00 

35.00 

35.00 

35.00 

35.00 

35.00 

35.00 

35.00 

35.00 

35.00 

35.00 

35.00 

35.00 

35.00 

35.00 

35.00 

35.00 

35.00 

35.00 

35.00 

35.00 


LONGITUDE 


0.00 

15.00 

30.00 

45.00 

60.00 

75.00 

90.00 

105.00 

120.00 

135.00 

150.00 

165.00 

180.00 

195.00 

210.00 

225.00 

240.00 

255.00 

270.00 

285.00 

300.00 

315.00 

330.00 

345.00 

0.00 

15.00 

30.00 

45.00 

60.00 

75.00 

90.00 

105.00 

120.00 

135.00 

150.00 

165.00 

180.00 

195.00 

210.00 

225.00 

240.00 

255.00 

270.00 

285.00 

300.00 

315.00 

330.00 

345.00 


RIGIDITY 

(GeV/ec) 

6.75 
7.27 

7.48 

7.70 

8.19 

8.73 
9.14 
9.29 

9.49 
9.89 

9.74 
8.95 
7.86 
6.46 

5.41 
4.55 
3.61 

2.76 
2.07 
1.93 

2.42 
3.41 
4.82 
5.92 
9.54 

9.89 
10.10 
10.53 
11.15 
11.44 
11.52 
11.71 
11.93 
12.04 
11.55 
10.60 

9.49 

8.97 
7.65 
6.12 
5.21 
4.25 

3.19 

2.89 
3.58 

4.98 
6.80 

8.70 


LATITUDE 


30.00 
30.00 
30.00 
30.00 
30.00 
30.00 
30.00 
30.00 
*  30.00 
30.00 
30.00 
30.00 
30.00 
30.00 
30.00 
30.00 
30.00 
30.00 
30.00 
30.00 
30.00 
30.00 
30.00 
30.00 
25.00 
25.00 
25.00 
25.00 
25.00 
25.00 
25.00 
25.00 
25.00 
25.00 
25.00 
25.00 
25.00 
25.00 
25.00 
25.00 
25.00 
25.00 
25.00 
25.00 
25.00 
25.00 
25.00 
25.00 


LONGITUDE 


0.00 

15.00 

30.00 

45.00 

60.00 

75.00 

90.00 

105.00 

120.00 

135.00 

150.00 

165.00 

180.00 

195.00 

210.00 

225.00 

240.00 

255.00 

270.00 

285.00 

300.00 

315.00 

330.00 

345.00 

0.00 

15.00 

30.00 

45.00 

60.00 

75.00 

90.00 

105.00 

120.00 

135.00 

150.00 

165.00 

180.00 

195.00 

210.00 

225.00 

240.00 

255.00 

270.00 

285.00 

300.00 

315.00 

330.00 

345.00 


RIGIDITY 

(GeV/eo) 

11.30 

11.71 
12.13 
12.67 
13.34 
14.07 
14.37 
14.40 
14.26 
13.95 

13.44 

12.72 
11,65 

10.48 
9.63 
8.78 
7.00 
5.60 

4.44 
4.07 
4.87 
6.98 
9.67 

10.60 

13.10 

13.64 

14.10 

14.53 

15.06 

15.58 

15.85 

15.79 

15.49 
15.03 

14.44 
13.76 
13.07 
12.43 
11.81 
10.98 

9.74 

7.89 

6.08 

5.44 
6.56 
9.65 

11.47 

12.46 
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LATITUDE 

LONGITUDE 

RIGIDITY 

LATITUDE 

LONGITUDE 

RIGIDITY 

(GeV/ec) 

(GeV/ec) 

20.00 

0.00 

14. 11 

10.00 

0.00 

14.73 

20.00 

15.00 

14.62 

10.00 

15.00 

15.26 

20.00 

30.00 

15.09 

10.00 

30.00 

15.80 

20.00 

45.00 

15.57 

10.00 

45.00 

16.34 

20.00 

60.00 

16.12 

10.00 

60.00 

16.94 

20,00 

75.00 

16.63 

10.00 

75.00 

17.44 

20.00 

90.00 

16.87 

10.00 

90.00 

17.67 

20.00 

105.00 

16.75 

10.00 

105.00 

17.56 

20.00 

120.00 

16.37 

*  10.00 

120.00 

17.18 

20.00 

135.00 

15.83 

10.00 

135.00 

16.65 

20.00 

150.00 

15.20 

10.00 

150.00 

16. 10 

20.00 

165.00 

14.55 

10.00 

165.00 

15.61 

20.00 

180.00 

13.93 

10.00 

180.00 

15.16 

20.00 

195.00 

13.40 

10.00 

195.00 

14.75 

20.00 

210,00 

12.91 

10.00 

210.00 

14.39 

20.00 

225.00 

12.33 

10.00 

225.00 

14.00 

20.00 

240.00 

11.47 

10.00 

240.00 

13.44 

20.00 

255.00 

9.69 

10.00 

255.00 

12.55 

■  20.00 

270.00 

7.82 

10.00 

270.00 

11.50 

20.00 

285.00 

6.84 

10.00 

285.00 

11.07 

20.00 

*  300.00 

8.65 

10.00 

300.00 

12. 16 

20.00 

315.00 

11.45 

10.00 

315.00 

13.09 

20.00 

330.00 

12.76 

10.00 

330.00 

'  13.76 

20.00 

345.00 

13.54 

10.00 

345.00 

14.25 

15.00 

0.00 

14.61 

5.00 

0,00 

14.50 

15.00 

15.00 

15.14 

5.00 

15.00 

14.99 

15.00 

30.00 

15.65 

5.00 

30.00 

15.52 

15.00 

45.00 

16.17 

5.00 

45.00 

16.10 

15.00 

60.00 

16.74 

5.00 

60.00 

16.71 

15.00 

75.00 

17.25 

5.00 

75.00 

17.22 

15.00 

90.00 

17.47 

5.00 

90.00 

17.47 

15.00 

105.00 

17.34 

5.00 

105.00 

17.42 

15.00 

120.00 

16.93 

5.00 

120.00 

17.11 

15.00 

135.00 

16.37 

5.00 

135.00 

16.66 

15.00 

150.00 

15.76 

5.00 

150.00 

16.21 

15.00 

165.00 

15.17 

5.00 

165.00 

15.83 

15.00 

180.00 

14.63 

5.00 

180.00 

15.49 

15.00 

195.00 

14.16 

5.00 

195.00 

15.15 

15.00 

210.00 

13.76 

5.00 

210.00 

14.82 

15.00 

225.00 

13.30 

5.00 

225.00 

14.47 

15.00 

240.00 

12.62 

5.00 

240.00 

14.03 

15.00 

255.00 

11.27 

5.00 

255.00 

13.41 

15.00 

270.00 

9.49 

5.00 

270.00 

12.67 

15.00 

285.00 

8.46 

5.00 

285.00 

12.40 

15.00 

300.00 

10.71 

5.00 

300.00 

12.78 

15.00 

315.00 

12.52 

5.00 

315.00 

13.36 

15.00 

330.00 

13.43 

5.00 

330.00 

13.80 

15.00 

345.00 

14.07 

5.00 

345.00 

14.11 
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LATITUDE 

LONGITUDE 

RIGIDITY 

(GeV/ec) 

LATITUDE 

LONGITUDE 

RIGIDITY 

(GeV/ec) 

0.00 

0.00 

13.94 

-10.00 

0.00 

12.11 

0.00 

15.00 

14.37 

-10.00 

15.00 

12.31 

0.00 

30.00 

14.87 

-10.00 

30.00 

12.71 

0.00 

45.00 

15.46 

-10.00 

45.00 

13.29 

0.00 

60.00 

16.10 

-10.00 

60.00 

13.89 

0.00 

75.00 

16.62 

-10.00 

75.00 

14.35 

0.00 

90.00 

16.90 

-10.00 

90.00 

14.68 

0.00 

105.00 

16.94 

-10.00 

105.00 

14.90 

0.00 

120.00 

16.73 

.  -10.00 

120.00 

14.91 

0.00 

135.00 

16.38 

-10.00 

135.00 

14.80 

0.00 

150.00 

16.05 

-10.00 

150.00 

14.74 

0.00 

165.00 

15.81 

-10.00 

165.00 

14.84 

0.00 

180.00 

15.59 

-10.00 

180.00 

14.94 

0.00 

195.00 

15.32 

-10.00 

195.00 

14.90 

0.00 

210.00 

15.03 

-10.00 

210.00 

14.76 

0.00 

225.00 

14.71 

-10.00 

225.00 

14.55 

0.00 

240.00 

14.33 

-10.00 

240.00 

14.29 

0.00 

255.00 

13.06 

-10.00 

255.00 

13.94 

0.00 

270.00 

13.32 

-10.00 

270.00 

13.51 

0.00 

285.00 

12.95 

-10.00 

285.00 

13.10 

0.00 

300.00 

13.05 

-10.00 

300.00 

12.92 

0.00 

315.00 

13.38 

-10.00 

315.00 

12.82 

0.00 

330.00 

13.58 

-10.00 

330.00 

12.55 

0.00 

345.00 

13.69 

-10.00 

345.00 

12.20 

-5.00 

0.00 

13.13 

-15.00 

0.00 

10.75 

-5.00 

15.00 

13.45 

-15.00 

15.00 

10.91 

-5.00 

30.00 

13.91 

-15.00 

30.00 

11.23 

-5.00 

45.00 

14.50 

-15.00 

45.00 

11.69 

-5.00 

60.00 

15.14 

-15.00 

60.00 

12.18 

-5.00 

75.00 

15.65 

-15.00 

75.00 

12.70 

-5.00 

90.00 

15.97 

-15.00 

90.00 

13.00 

-5.00 

105.00 

16. 10 

-15.00 

105.00 

13.25 

-5.00 

120.00 

16.00 

-15.00 

120.00 

13.25 

-5.00 

135.00 

15.77 

-15.00 

135.00 

13.32 

-5.00 

150.00 

15.58 

-15.00 

150.00 

13.48 

-5.00 

165.00 

15.50 

-15.00 

165.00 

13.78 

-5.00 

180.00 

15.42 

-15.00 

180.00 

14.09 

-5.00 

195.00 

15.25 

-15.00 

195.00 

14.24 

-5.00 

210.00 

15.02 

-15.00 

210.00 

14.24 

-5.00 

225.00 

14.74 

-15.00 

225.00 

14. 15 

-5.00 

240.00 

14.41 

-15.00 

240.00 

13.97 

-5.00 

255.00 

14.01 

-15.00 

255.00 

13.69 

-5.00 

270.00 

13.53 

-15.00 

270.00 

13.30 

-5.00 

285.00 

13.14 

-15.00 

285.00 

12.89 

-5.00 

300.00 

13.08 

-15.00 

300.00 

12.59 

-5.00 

315.00 

13.19 

-15.00 

315.00 

12.31 

-5.00 

330.00 

13.15 

-15.00 

330.00 

11.76 

-5.00 

345.00 

13.04 

-15.00 

345.00 

11.15 
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LATITUDE 

LONGITUDE 

RIGIDITY 

(GeV/ec) 

-20,00 

0.00 

9.21 

-20.00 

15.00 

9.06 

-20,00 

30.00 

9.29 

-20.00 

45.00 

9.70 

-20,00 

60.00 

10.21 

-20.00 

75.00 

10.25 

-20.00 

90.00 

10.41 

-20.00 

105.00 

10.65 

-20.00 

120.00 

10.84 

-20.00 

135.00 

10.62 

-20.  00 

150.00 

10.74 

-20.00 

165.00 

11.93 

-20.00 

180.00 

12.78 

-20.00 

195.00 

13.19 

-20.00 

210.00 

13.44 

-20.00 

225.00 

13.51 

-20.00 

240.00 

13.46 

-20.00 

255.00 

13.28 

-20.00 

270.00 

12.94 

-20.00 

285.00 

12.53 

-20.00 

300.00 

12.14 

-20.00 

315.00 

11.65 

-20.00 

330.00 

*10.73 

-20.00 

345.00 

9.78 

-25.00 

0.00 

7.55 

-25.00 

15.00 

7.50 

-25.00 

30.00 

7.72 

-25.00 

45.00 

7.88 

-25.00 

60.00 

7.86 

-25.00 

75.00 

7.37 

-25.00 

90.00 

7.08 

-25.00 

105.00 

7.24 

-25.00 

120.00 

7.43 

-25.00 

135.00 

7.71 

-25.00 

150.00 

8.47 

-25.00 

165.00 

9.55 

-25.00 

180.00 

10.22 

-25.00 

195.00 

11.22 

-25.00 

210.00 

12.01 

-25.00 

225.00 

12.63 

-25.00 

240.00 

12.76 

-25.00 

255.00 

12.72 

-25.00 

270.00 

12.46 

-25.00 

285.00 

12.04 

-25.00 

300.00 

11.53 

-25.00 

315.00 

10.74 

-25.00 

330.00 

9.63 

-25.00 

345.00 

8.36 

LATITUDE 

LONGITUDE 

RIGIDITY 

(GeV/ec) 

-30.00 

0.00 

6.36 

-30.00 

15.00 

6.02 

-30.00 

30.00 

5.86 

-30.00 

45.00 

5.79 

-30.00 

60.00 

5.42 

-30.00 

75.00 

5.23 

-30. 00 

90.00 

5.10 

-30.00 

105.00 

5.16 

-30.00 

120.00 

5.19 

-30.00 

135.00 

5.39 

-30.00 

150.00 

5.89 

-30.00 

165.00 

6.58 

-30.00 

180.00 

7. 99 

-30.00 

195.00 

9.45 

-30.00 

210.00 

9,  43 

-30.00 

225.00 

10.64 

-30.00 

240.00 

11.85 

-30.00 

255.00 

12.00 

-30.00 

270.00 

11.88 

-30.00 

285.00 

11.43 

-30.00 

300.00 

10.75 

-30.00 

315.00 

9.83 

-30.00 

330.00 

8. 39 

-30.00 

345.00 

7.09 

-35.00 

0.00 

5.24 

-35.00 

15.00 

4.59 

-35.00 

30.00 

4.45 

-35.00 

45.00 

4.31 

-35.00 

60.00 

4.07 

-35.00 

75.00 

3.72 

-35.00 

90.00 

3.35 

-35.00 

105.00 

3.32 

-35.00 

120.00 

3.37 

-35.00 

135.00 

3.65 

-35.00 

150.00 

4. 10 

-35.00 

165.00 

4.90 

-35.00 

180.00 

5.65 

-35.00 

195.00 

6.54 

-35.00 

210.00 

7.88 

-35.00 

225.00 

9.11 

-35.00 

240.00 

9.58 

-35.00 

255.00 

11.12 

-35.00 

270.00 

11.16 

-35.00 

285.00 

10.67 

-35.00 

300.00 

9.90 

-35.00 

315.00 

8.72 

-35.00 

330.00 

7. 16 

-35.00 

345.00 

6.18 
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LATITUDE 

LONGITUDE 

RIGIDITY 

LATITUDE 

LONGITUDE 

RIGIDITY 

(GeV/ec) 

(GeV/ec) 

-40.00 

0.00 

4.29 

-50.00 

0.00 

2.87 

-40.00 

15.00 

3.74 

-50.00 

15.00 

2.37 

-40.00 

30.00 

3.40 

-50.00 

30.00 

1.95 

-40.00. 

45.00 

3.21 

-50.00 

45.00 

1.60 

-40.00 

60.00 

2. 83 

-50.00 

60.00 

1.27 

-40.00 

75.00 

2.43 

-50.00 

75.00 

0.94 

-40.00 

90.00 

2.08 

-50.00 

90.00 

0.66 

-40.00 

105.00 

2.00 

-50.00 

105.00 

0.54 

-40.00 

120.00 

2.01 

.  -50.00 

120.00 

0.53 

-40.00 

135.00 

2.21 

-50.00 

135.00 

0.60 

-40.00 

150.00 

2.65 

-50.00 

150.00 

0.84 

-40.00 

165.00 

3.24 

-50.00 

165.00 

1.15 

-40.00 

180.00 

4.11 

-50.00 

180.00 

1.63 

-40.00 

195.00 

4.76 

-50.00 

195.00 

2.24 

-40.00 

210.00 

5.56 

-50.00 

210.00 

2.94 

-40.00 

225.00 

6.65 

-50.00 

225.00 

3.76 

-40.00 

240.00 

8.20 

-50.00 

240.00 

4.48 

-40.00 

255.00 

9.75 

-50.00 

255.00 

5.57 

-40.00 

270.00 

10.18 

-50.00 

270.00 

7.02 

-40.00 

285.00 

9.77 

-50.00 

285.00 

7.57 

-40.00 

300.00 

8.98 

-50.00 

300.00 

6.98 

-40.00 

315.00 

7.61 

-50.00 

315.00 

5.68 

-40.00 

330.00 

6.42 

-50.00 

330.00 

4.51 

-40.00 

345.00 

5.31 

-50. 00 

345.00 

3.50 

-45.00 

0.00 

3.46 

-55.00 

0.00 

2.23 

-45.00 

15.00 

2.96 

-55.00 

15.00 

1.79 

-45.00 

30.00 

2.53 

-55.00 

30.00 

1.42 

-45.00 

45.00 

2.30 

-55.00 

45.00 

1.12 

-45.00 

60.00 

1.93 

-55.00 

60.00 

0.87 

-45.00 

75.00 

1.53 

-55.00 

75.00 

0.53 

-45.00 

90.00 

1.28 

-55.00 

90.00 

0.34 

-45.00 

105.00 

1.12 

-55.00 

105.00 

0.23 

-45.00 

120.00 

1.11 

-55.00 

120.00 

0.21 

-45.00 

135.00 

1.25 

-55.00 

135.00 

0.26 

-45.00 

150.00 

1.51 

-55.00 

150.00 

0.38 

-45.00 

165.00 

2.04 

-55.00 

165.00 

0.59 

-45.00 

180.00 

2.72 

-55.00 

180.00 

0.90 

-45.00 

195.00 

3.33 

-55.00 

195.00 

1.38 

-45.00 

210.00 

4.24 

-55.00 

210.00 

1.88 

-45.00 

225.00 

4.93 

-55.00 

225.00 

2.64 

-45.00 

240.00 

5.91 

-55.00 

240.00 

3.38 

-45.00 

255.00 

7.83 

-55.00 

255.00 

4.20 

-45.00 

270.00 

9.00 

-55.00 

270.00 

4.96 

-45.00 

285.00 

8.76 

-55.00 

285.00 

5. 19 

-45.00 

300.00 

7.84 

-55.00 

300.00 

5.02 

-45.00 

315.00 

6.91 

-55.00 

315.00 

4.45 

-45.00 

330.00 

5.63 

-55.00 

330.00 

3.67 

-45.00 

3^5.00 

4.32 

-55.00 

345.00 

2.93 
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LATITUDE 

LONGITUDE 

RIGIDITY 

(GeV/ec) 

-60.00 

0.00 

1.78 

-60.00 

15.00 

1.32 

-60.00 

30.00 

1.03 

-60.00 

45.00 

0.75 

-60.00 

60.00 

0.49 

-60.00 

75.00 

0.30 

-60.00 

90.00 

0.15 

-60.00 

105.00 

0.08 

-60.00 

120.00 

0.06 

-60.00 

135.00 

0.08 

-60.00 

150.00 

0. 14 

-60.00 

165.00 

0.27 

-60.00 

180.00 

0.48 

-60.00 

195.00 

0.79 

-60.00 

210.00 

1.18 

-60.00 

225.00 

1.62 

-60.00 

240.00 

2.23 

-60.00 

255.00 

3.00 

-60.00 

270.00 

3.77 

-60.00 

285.00 

3.95 

-60.00 

300.00 

3.97 

-60.00 

315.00 

3.52 

-60.00 

330.00 

2.88 

-60.00 

345.00 

2,27 

-65.00 

0.00 

1.30 

-65.00 

15.00 

0.98 

-65.00 

30.00 

0.72 

-65.00 

45.00 

0.50 

-65.00 

60.00 

0.30 

-65.00 

75.00 

0.15 

-65.00 

90.00 

0.06 

-65.00 

105.00 

0.00 

-65.00 

120.00 

0.00 

-65.00 

135.00 

0.00 

-65.00 

150.00 

0.03 

-65.00 

165.00 

0.11 

-65. 00 

180.00 

0.23 

-65.00 

195.00 

0.43 

-65.00 

210.00 

0.71 

-65.00 

225.00 

1.06 

-65.00 

240.00 

1.51 

-65.00 

255.00 

1.98 

-65.00 

270.00 

2.53 

-65.00 

285.00 

2.71 

-65.00 

300.00 

2.72 

-65.00 

315.00 

2.50 

-65.00 

330.00 

2. 10 

-65.00 

345.00 

1.61 

LATITUDE 

LONGITUDE 

RIGIDITY 

(GeV/ec) 

-70.00 

0.00 

0.89 

-70.00 

15.00 

0.64 

-70.00 

30.00 

0.47 

-70.00 

45.00 

0.31 

-70.00 

60,00 

0.18 

-70.00 

”5.00 

0.08 

-70.00 

90.  t.  0 

0.00 

-70.00 

ior.n; 

0.00 

-70.  JO 

'  i.u,  1 

0.00 

-70  ,0 

135,00 

0.00 

-70.00 

1‘>'\  00 

0.00 

-70.00 

165.00 

0  J3 

-70.00 

180.00 

c,  10 

-70.00 

1J0.00 

0.22 

-70.00 

210.00 

0.41 

-70.70 

225  00 

0.S4 

-70.00 

240.00 

0.  96 

-70.00 

255. 00 

-70. 00 

270.00 

1.58 

-70.00 

285.00 

1,75 

-70.00 

300.00 

1.80 

-70.00 

315.00 

1.67 

-70.00 

330.00 

1.39 

-70.00 

345.00 

1.14 

-75.00 

0.00 

0.59 

-75.00 

15.00 

0.43 

-75.00 

30.00 

0.30 

-75.00 

45.00 

0.19 

-75.00 

60.00 

0.10 

-75.00 

75.00 

0.04 

-75.00 

90.00 

0.00 

-75.00 

105.00 

0.00 

-75.00 

120.00 

0.00 

-75.00 

135.00 

0.00 

-75.00 

150.00 

0.00 

-75.00 

165.00 

0.00 

-75. 00 

180.00 

0.05 

-75. 00 

195.00 

0.12 

-75,00 

210.00 

0.23 

-75.00 

225.00 

0.36 

-75.00 

240.00 

0.54 

-75.00 

255.00 

0.72 

-75.00 

270.00 

0.91 

-75.00 

285.00 

1.02 

-75.00 

300.00 

1.03 

-75.00 

315.00 

1.05 

-75.00 

330.00 

0.88 

-75.00 

345.00 

0.76 
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EM 


LATITUDE 

LONGITUDE 

RIGIDITY 

(GeV/ec) 

-80.00 

0.00 

0.37 

-80.00 

15.00 

0.28 

-80.00 

30.00 

0.19 

-80.00 

45.00 

0. 13 

-80.00 

60.00 

0.07 

-80.00 

75.00  . 

0.03 

-80.00 

90.00 

0.00 

-80.00 

105.00 

0.00 

-80.00 

120.00 

0.00 

-80.00 

135.00 

0.00 

-80.00 

150.00 

0.00 

-80.00 

165.00 

0.00 

-80.00 

180.00 

0.03 

-80.00 

195.00 

0.08 

-80.00 

210.00 

0.13 

-80.00 

225.00 

0.21 

-80.00 

240.00 

0.30 

-80.00 

255.00 

0.40 

-80.00 

270.00 

0.48 

-80.00 

285.00 

0.55 

-80.00 

300,00 

0.56 

-80.00 

315.00 

0.54 

-80.00 

330.00 

0.49 

-80.00 

345.00 

0.42 

oooo  ooonooooooooooooooo 


APPENDIX  7:  SOURCE  CODES  OF  THE  CREME  PROGRAMS 


PROGRAM  STASS 

THIS  PROGRAM  PREPARES  ORBIT -AVERAGED  TRAPPED  PROTON  DATA 
FILES  FOR  THE  LET  PROGRAM. 

THE  TRAPPED  PROTON  DATA  CAN  BE  OBTAINED 
FROM  REPORTS  PUBLISHED  BY  THE  NATIONAL  SPACE  SCIENCE  DATA 
CENTER  AT  GODDARD  SPACE  FLIGHT  CENTER.  THE 
PROGRAM  IS  SET  UP  TO  TAKE  THE  TABULAR  DATA  UNDER 
'AVERAGED  DIFFER.  FLUX'  IN  THE  'COMPOSITE  ORBIT 
SPECTRUM'  TABLE  OF  ONE  OF  THESE  REPORTS.  CARE  SHOULD  BE 
TAKEN  TO  BE  SURE  CONTRIBUTIONS  FROM  SOLAR  FLARE  PROTONS 
HAVE  NOT  BEEN  INCLUDED  IN  THIS  TABLE. 

THE  DATA  SHOULD  BF  ENTERED  IN  PAIRS  OF  POINTS  (I.E.  THE 
PROTON  ENERGY  FOLLOWED  BY  THE  CORRESPONDING  FLUX). 

THE  ENERGY  SHOULD  BE  IN  MEV  AND  THE  FLUX  SHOULD  BE  IN 
#/CM**2/SEC/KEV.  DATA  BELOW  10  MEV  NEED  NOT  BE  ENTERED. 

ENTER  A  PAIR  OF  ZERO  VALUES  TO  CLOSE  THE  FILE  AFTER  ALL 
THE  DATA  HAS  BEEN  ENTERED.  THIS  FILE  MAY  BE  EDITED 
AFTERWARDS  TO  CORRECT  ANY  ERRORS  IN  THE  DATA  YOU  ENTERED. 

CHARACTERS  ALPHA 
WRITE  (6, 100) 

100  FORMAT ('  ENTER  THE  NAME  OF  THE  FILE  TO  CONTAIN  THE  PROTON  DATA 

1:  ',$> 

READ  (5  , 200  )ALPHA 

NOTE?  THIS  FILE  MUST  BE  RENAMED  TO  STASS. DAT  BEFORE  IT 
CAN  BE  USED  IN  THE  LET  PROGRAM. 

200  FORMAT (A 12) 

OPEN (UNIT=1 ,  STATUS  =  'NEW'  ,FILE=ALPHA ) 

WRITE  (6, 300) 

300  FORMAT ('  ENTER  THE  DATA  POINT  PAIRS,  STARTING  WITH '/ 

2  '  10  MEV  AND  PROCEEDING  MONOTONICALLY  TO  THE  HIGHEST'/ 

4  '  ENERGY.  THE  PAIRS  ARE  TO  BE  TAKEN  FROM  THE  COMPOSITE'/ 

5  '  ORBIT  SPECTRUM.  THE  PAIRS  ARE  ENERGY  LEVEL  IN  MEV,  FOLLOWED'/ 

6  '  BY  THE  AVERAGED  DIFFER.  FLUX  IN  0/CM»«2/SEC/KEV. ' / 

7  '  ENTER  A  PAIR  OF  ZERO  VALUES  WHEN  YOU  ARE  DONE.') 

1  =  1 

6  CONTINUE 
WRITE  (6,500)1 

500  FORMAT ('  ENTER  PAIR  ',12, »:  ',$) 

ACCEPT  * ,E,F 
IF  (E.LE.O.)  GOTO  7 
WRITE  (1,600  )E,F 

600  FORMAT  (1X,F6. 1, 1X,E10.3) 

1=1+1 
GOTO  6 

7  CONTINUE 
END 


ooooooo  o  o  o  o  o  onooo  OOOOOOOOOOOOOOOOOOO 


PROGRAM  GEOMAG 

THIS  PROGRAM  CALCULATES  THE  2  DAY  ORBIT  AVERAGE  OF 

THE  GEOMAGNETIC  TRANSMISSION  FUNCTION  FOR  CIRCULAR  ORBITS. 

WHEN  YOU  RUN  THIS  PROGRAM,  ALL  INPUT 

DATA  WILL  BE  REQUESTED  IN  PROMPTS.  OUTPUT  DATA  TABULATED  IN 
INTERVALS  OF  . 1  GV  ARE  STORED  IN  THE  OUTPUT  FILE  GTRANS.DAT. 
[EXCEPTION:  THE  HIGHEST  RIGIDITY  VALUE  WITH  A  ZERO 
TRANSMISSION  FUNCTION  (I.E. ,  A  PARTICLE  WITH  THAT  RIGITITY 
CANNOT  REACH  THE  SPACECRAFT  AT  ANY  POINT  ON  THE  SPACECRAFT’S 
ORBIT)  MAY  BE  INTERPOLATED  TO  GIVE  A  VALUE  NOT  DIVISIBLE  BY 
. 1  GV.  THE  RIGIDITIES  WILL  BE  IN  INCREASING  ORDER  IN  ANY  CASE. ] 
AN  INPUT  FILE  CALLED  " CUTOFF. DAT"  CONTAINS  THE  TABULATION 
OF  WORLDWIDE  VERTICAL  GEOMAGNETIC  CUTOFFS  AT  20KM  ALTITUDE, 

TAKEN  FROM  M.  A.  SHEA  AND  D.  F.  SMART,  REPORT  NO. 

AFCRL-TR -75-01 85,  HANSCOM  AFB,  MASS.,  1975. 

AZ  IS  THE  AZMUTH  ANGLE  OF  THE  PARTICLE  WRT  THE  SPACECRAFT. 

ZE  IS  THE  ZENITH  ANGLE  OF  THE  PARTICLE  WRT  THE  SPACECRAFT. 

THESE  ANGLES  ARE  BOTH  IN  DEGREES  AND  ARE  DATAED  IN  AS  ZERO. 

THEY  CAN  BE  RESET  TO  ANY  OTHER  ARRIVAL  DIRECTION. 

DIMENSION  MAT (200 )  ,CUTOFF(33, 24 ) ,T(201 ) ,CF(201 ) 

DATA  ( MAT  (J ) ,  J  =  1 , 200 )  /200  *0  / 

DATA  AZ/O .  /  ,ZE/0 .  /  ,AZG/0 .  /  ,ZEG/0 .  / 

OPEN (UNIT=1 5, READONLY, SHARED, STATUS* ’OLD’, FILE= ’CUTOFF. DAT’) 

OPE N  (UNIT =1  6 ,  FILE  a  'GTR ANS .  DAT  ’ ,  STATUS  a  'NEW  ’  ) 

ASK  FOR  INPUT. 

IS  THE  SHADOW  OF  THE  EARTH  ON  THE  SPACECRAFT 
TO  BE  ACCOUNTED  FOR  IN  THIS  CALCULATION? 

WRITE  (6,411) 

ACCEPT  «,  ISHADOW 

IS  THE  MAGNETOSPHERE  QUIET  OR  STORMY? 

WRITE  (6,412) 

ACCEPT  *,  ISTORM 


CALL  ORBIT  (1, PERIOD, ZLON,ZLAT, RADIUS) 

WRITE  (6,1)  PERIOD 

FORMAT (’  THE  ORBITAL  PERIOD  IS  ',E12.5,'  SECONDS.’) 

READ  IN  THE  TABLE  OF  WORLD  WIDE  VERTICAL 
GEOMAGNETIC  CUTOFFS  AT  20KM  ALTITUDE. 

THESE  DATA  ARE  TABULATED  EVERY  5  DEGREES  IN  LATITUDE 
AND  EVERY  15  DEGREES  IN  LONGITUDE. 

THEY  ARE  TAKEN  FROM  M.  A.  SHEA  AND  D.  F.  SMART,  REPORT  NO.  > 

AFCRL-TR -75-01 85,  HANSCOM  AFB,  MASS,  1975. 

DO  20  I  si ,  33 
11=34-1 


o  o  o  o  ooo  oooo  o  o  o  oooo  oooo  o  o  o 


DO  20  Js1,24 

READ<15,405)  CUTOFFdl.J)  /  . 

THIS  ACCOUNTS  FOR  GEOMAGNElIC  pUTOFF  SUPRESSION  DURING 
LARGE  MAGNETIC  STORMS  (FOLLOWING  ADAMS  ET  AL. ,  1981). 

IF(ISTORM.EQ.  1 )  CUTOFFdl, J|MC:^::>jFdI, J)*(1  .-.54»ZXP< 

1  -CUT0FF(IIf J)/2,9))  V!'  • 

CONTINUE  V'/ 

AVERAGE  THE  CUTOFF  AROUND  THE  80  DEGREE  LATITUDE  LINES, 
NORTH  AND  SOUTH,  TO  APPROXIMATE  THE  CUTOFFS  AT  90  DEGREES. 


CN=0 . 

CS=0. 

DO  30  J  si ,  24 
CN=CN*CUTOFF(33,J) 

30  CS=CS4CUT0FF(1,J) 

CNsCN/24. 

CSsCS/24. 

COMPUTE  THE  TOTAL  NUMBER  OF  STEPS  IN  TWO  DAYS  IF  WE  MAKE 
200  STEPS  PER  ORBIT. 

JMAXsINT (2. *200. *86400. /PERIOD  +1.5) 

COMPUTE  THE  STEP  SIZE  IN  SECONDS. 

STEP=PERIOD/200. 

COMPUTE  THE  VERTICAL  CUTOFF  AT  THE  SPACECRAFT 
POSITION  FOR  EVERY  TIME  STEP. 

DO  50  J=1 ,  JMAX 
TIME  sFLOAT (J -1 ) *STEP 
CALL  ORBIT (2, TIME, ZL0N,ZLAT, ALT) 

COMPUTE  THE  TABULAR  POSITION  OF  THE  VERTICAL  CUTOFF. 

ZIsZLAT/5.  +17. 

ZJ=ZLON/15.  +  1. 

ILOsINT(ZI) 

IUPsILO+1 
JLOsINT(ZJ) 

JUPsJLO+1 

IF  (JUP.  EQ.  25)  JUPsI 

INTERPOLATE  THE  VERTICAL  CUTOFF  TO  THE  EXACT  LOCATION 
OF  THE  SPACECRAFT  USING  STORMER  THEORY. 

IF( ABS (ZLAT ) .GE. 80. )  GO  TO  100 
DIsZI -FLOAT  (ILO) 

DJsZJ-FLOAT(JLO) 
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XORBs(6371 . 2+ALT  )/6371 . 2 
RGRD*  6391.2/6371.2 
ZL0NLQs(JL0-1  )*1 5. 

ZLONUP  =Z  LONLO+1 5.  • 

ZLATLO=ILO«5.-85. 

Z  LATUP  aZLATLO+5, 

SC  sSTORMER ( ZLAT, ZLON , XORB , AZ , ZE ) 

SCLL sSTORMER ( Z  LATLO , Z  LONLO ,  R  GR  D,  AZ  G,  ZEG ) 

SC UL sSTORMER (ZLATUP, ZLONLO, RGRD, AZG, ZEG ) 

SCLU sSTORMER (Z LATLO,  ZLONUP,  RGRD,  AZG,  ZEG)  . 

SCUU sSTORMER (ZLATUP,  ZLONUP,  RGRD,  AZG,  ZEG) 

Y1 =SC*CUTOFF(ILO, JLO)/SCLL 
Y2s  SC# CUTOFF (I UP, JLO )/SCUL 
Y3s  SC#C  UTOFF ( ILO , JUP ) / SC  LU 
Y4s  SC* CUTOFF (I UP, JUP)/SCUU  " 

Cs  ( 1 .  ~DI )  *  ( 1 .  -DJ  )  *Y1  +  ( 1 .  -DI )  *DJ  *Y3  +DI  *  ( 1 .  -DJ )  *Y2+DI  *DJ  *Y  4 
GO  TO  200 
C 

C  FOR  ABS (LATITUDE  ).GT. 80  USE  THE  CUTOFFS  AT  THE  POLE  INSTEAD  OF 

C  THE  CUTOFFS  AT  FOUR  NEARBY  LOCATIONS. 

C 

100  CONTINUE 

IF( ZLAT.LE. -80. )G0  TO  110 
SC  sSTORMER (90.,0.,X0RB,AZ,ZE> 

SCXsSTORMER  (90.  ,0. , RGRD,  AZG,  ZEG) 

C=SC*CN/SCX 
GO  TO  200 
110  CONTINUE 

SC sSTORMER (-90., 0. ,XORB,AZ,ZE) 

SCX  sSTORMER (-90. , 0. , RGRD, AZG, ZEG ) 

CsSC*CS/SCX 
200  CONTINUE 

C 

C  HISTOGRAM  THE  CUTOFFS. 

C 

IDEXsINT (C*1 0. )+1 
MAT (IDEX )=MAT (IDEX )+1 

C  TO  PRINT  OUT  THE  LONGITUDE,  LATITUDE,  AND  CUTOFF  AT  EACH 

C  POINT  ARROUND  THE  SPACECRAFT’S  ORBIT,  JUST  ENABLE  THE  FOLLOWING 

C  STATEMENT  BY  REMOVING  THE  COMMENT  C. 

C  WRITE (6, 435 )ZLON, ZLAT, C 

50  CONTINUE 

CMAT=0. 

C 

C  SAVE  THE  THRESHOLD. 

C 

JSAVsO 

DO  300  J si ,200 
C 

C  CONVERT  THE  HISTOGRAM  TO  TRANSMISSION. 


o  o  o  o  o  o  o 


CMATsFLOAT (MAT (J  ) ) /FLOAT (JHAX )+CMAT 
IF(  (JSAV.EQ.O).AND.(CMAT.GT.O))  JSAVsJ+1 
CFCJ+1  )«FLOAT(J  )/10. 

300  TU+1 )  =CMAT 
CF( 1 )*0, 0 

T  (1  )«0.0 

INTERPOLATE  THE  HIGHEST  ZERO  TRANSMISSION  VALUE  TO  BE 
JUST  AT  THE  CUTOFF  THRESHOLD. 

DTsT  (JSAV+1  )-T  (JSAV) 

IF(DT.EQ. 0. )  GO  TO  301 

CFTsCF(JSAV>-T  (JSAV)*(CF( JSAV+1  )-CF(JSAV)  )/DT 
CF(JSAV-1  )=AMAX1(CFT,CF(JSAV-1  )) 

301  CONTINUE 

THIS  IS  A  CORRECTION  FOR  THE  EARTH’S  SHADOW  ON  THE  SPACECRAFT 
ACCORDING  TO  SIMPLE  GEOMETRICAL  OPTICS. 

IF(ISHADOW.NE.I)  GO  TO  390 
DO  380  J  =  1,200 

T (J  )=T  (J)«(1.-0.5»a.-<(6371.24ALT)«*2.-(6371.2)**2.)**.5/ 

1  (6371.2+ALT))) 

380  CONTINUE 

390  DO  400  J*1,200 

400  WRITE  (1  6,  410  )CF(  J  )  #T(J  ) 

405  FORMAT  (17X.F6.2) 

410  FORMAT (5X,  F6.  3,  5X,F8. 6) 

411  FORMAT (IX, 'DO  YOU  WANT  TO  INCLUDE  THE  EFFECT  OF  THE  SHADOW',/, 

1  '  OF  THE  EARTH?  (1  FOR  YES;  0  FOR  NO):  ',$) 

412  FORMAT (IX, 'ENTER  THE  MAGNETIC  WEATHER  CONTITION:  1  FOR  STORMY; 
1  0  FOR  QUIET:  ',$) 

435  FORMAT ('  LONGITUDE  =  »,F6.2,',  LATITUDE  =  ',F5.2,'  CUTOFF  *  ', 

1  F6.  3,'  . ' ) 

END 
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FUNCTION  STORMER (GCLATD,GCLOND,RGC,AZ,ZE ) 

C  WE  DID  NOT  WRITE  THIS  SUBROUTINE.  WE  HAVE  MADE  NO  CHANGES  IN 

C  IT  IN  1984.  C 

COMMON /KARL/RED, EDLAT, AZM, ZEM .GAMMA 

C  THIS  FUNCTION  TRANSFORMS  A  GEOGRAPHIC  LOCATION  AND  ARRIVAL 

C  DIRECTION  INTO  OFFSET  DIPOLE  COORDINATES,  THEN  COMPUTES  THE 

C  STORMER  CUTOFF  IN  GV  AND  RETURNS  THE  RESULT.  THE  OFFSET  DIPOLE 

C  COORDINATES  ARE  AVAILABLE  IN  THE  COMMON  BLOCK  /KARL/. 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


c 

c 

c 


c 


GCLATD  IS  GEOCENTRIC  LATITUDE  IN  DEGREES 

GCLONG  IS  GEOCENTRIC  LONGITUDE  IN  DEGREES 

RGC  IS  RADIAL  DISTANCE  FROM 'GEOCENTER  IN  EARTH  RADII 

AZ  IS  GEOGRAPHIC  AZIMUTH 

ZE  IS  GEOGRAPHIC  ZENITH 

RED  IS  RADIAL  DISTANCE  FROM  OFFSET  DIPOLE  POSITION  IN 

EARTH  RADII 

EDLAT  IS  THE  GEOMAGNETIC  LATITUDE  IN  OFFSET  DIPOLE  COORDINATES 
AZM  IS  GEOMAGNETIC  AZIMUTH  IN  OFFSET  DIPOLE  COORDINATES 

ZEM  IS  GEOMAGNETIC  ZENITH  IN  OFFSET  DIPOLE  COORDINATES 

GAMMA  IS  GAMMA  ANGLE  MEASURED  FROM  MAGNETIC  EAST 


DATA  JDATA ,  NOPT/2 *0  / , PI ,  RAD,  PI02,  XEDFGC ,  YEDFGC,  ZEDFGC , CP, SP 
1  ,ST,CPCT,CPST,SPCT,SPST,XGMED,YGMED,ZGMED/16*-8000./ 

DATA  SM ALL/1.  OE -35/ 

DPEC(U  )»SIGN  (1 . /SNGL  (DSQRTODO+DBLE  ( (U  /ZEDRTL)**2  ) ) )  ,ZRTL*ZEDRTL) 
IF( JDATA. EQ. 77 )  GO  TO  10 


PI  =  ACOS(-I.O) 

RAD  a  180.0 /PI 
PI02  a  PI/2.0 
TWOPI  a  PI  *2. 0 

DATA  ERAD,  THETAD,  PHID,  R1KM  ,TH  1DEG,  PH  1DEG/6371 . 2, 
1  11.4354,-290.2392,450.2586,72.8278,148.7753/ 
NOPT  a  0 

SQRT3  a  SORT  (3.0) 

ENTER  GEOMAGNETIC  DATA,  IGRF  1975 


DT  a  5 

SEE  JGR, 

.0 

81, 

5163,  1976 

G01  a 

-30186.0 

+ 

25. 6#DT 

G02  a 

-1  898.0 

- 

24.9*DT 

Gil  a 

-2036.0 

+ 

10.  0*DT 

G12  a 

2997.0 

+ 

0.7*DT 

H11  a 

5735.0 

- 

10. 2*DT 

H  12  S 

-2124.0 

3.0*DT 

G22  a 

1551.0 

+ 

4 . 3*DT 

H22  a 

-37.0 

- 

1  8.9*DT 

DT  IS  NUMBER  OF  YEARS  SINCE  1975 


IFCNOFT. EQ. 1 )PRINT  1000,  G01,  G02,  G11,  G12,  G22,  H1 1 ,  H12,  H22 
COMPUTE  POSITION  OF  OFFSET  DIPOLE 
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HO  a  SQRT  (  GOI  *G01  -+-Q1 1  *G1 1 +H 1 1  «H  1 1  ) 

HOSQ  a  HO*HO 

ELO  a  2. 0*G01  *G02+(G  1 1  *G  12+H 1 1  *H  12  )*SQRT3 
ELI  a  -G1 1  *G02+  CGOI  *G12+G11  *G22+H1 1  *H22 )  *SQRT3 
EL2  a  -H 1 1  *G02+(G01  *H  12-H 1 1  *G22-KJ  1 1  *H22)*SQRT3 
C  E  a  (EL0*G01  +EL1  *G1 1 +EL2  *H1 1 )  *4.0  *HOSQ 

E  a  (EL0*G01  +EL1  *G  1 1+EL2*H  1 1  )/(4. 0*HQSQ) 

IF(NOPT.EQ.I)  PRINT  1011,  ELO,  ELI,  EL2,  E,  HO 
1011  FORMAT (1H  ,  8E15.5) 

C  XEDFGC  a  ERAD*(EL1  -G1 1  *E  )/(3.0  *HOSQ ) 

XEDFGC  =  <EL1-G11*E)/(3.0*H0SQ) 

C  YEDFGC  a  ERAD*(EL2-H11*E)/(3.0«H0SQ) 

YEDFGC  a  (EL2-H11*E)/(3.0«H0SQ) 

C  ZEDFGC  a  ERAD*  (EL0-G01  *E)/(3.0*H0SQ ) 

ZEDFGC  a  (EL0-G01  *E  )/(3.  Q*H0SQ) 

REDFGC  a  SQRT(XEDFGC*XEDFGC+YEDFGC*YEDFGC+ZEDFGC*ZEDFGC) 
IF(NOPT.EQ.  DPRINT  3001,  XEDFGC,  YEDFGC,  ZEDFGC,  REDFGC 

3001  FORMAT  (1H  ,  4F10.4,  3X,  'XEDFGC,  YEDFGC,  ZEDFGC,  REDFGC') 

1000  FORMAT  (1H  ,  1 0F  1 3 « 5 ) 

1010  FORMAT  (1H0,  8F15.5/1H  ,  8F15.5  > 

THETA  a  THETAD/RAD 
PHI  a  PHID/RAD 
CP  a  COS  (PHI) 

SP  a  SIN  (PHI) 

ST  a  SIN  (THETA) 

CT  a  COS  (THETA) 

CPCT  a  CP*CT 
CPST  a  CP*ST 
SPCT  a  SP*CT 
SPST  a  SP *ST 
R1ER  a  R1KM/ERAD 
THI RAD  a  TH1DEG/RAD 
PH1RAD  a  PH  1DEG/RAD 

IF(NOPT.EQ.  1  )PRINT  1000,  R1KM ,TH1DEG,  PH  1DEG,  R1ER.TH  1RAD,  PH  1RAD 

XGMED  a  XEDFGC* CPCT  -YEDFGC*SPCT  -ZEDFGC*ST 

YGMED  a  XEDFGC*SP  +YEDFGC*CP 

ZGMED  a  XEDFGC*CPST  -YEDFGC*SPST  +ZEDFGC*CT 

IF(NOPT.EQ.  DPRINT  3002,  XGMED,  YGMED,  ZGMED 

3002  FORMAT (1H  ,  3F10.4,  13X,  'XGMED,  YGMED,  ZGMED') 

IF(NOPT.EQ. DPRINT  1010,  CP,  SP,  CT,  ST,  CPCT,  CPST,  SPCT,  SPST 
JDATA  a  77 

10  CONTINUE 

ITERATE  TO  FIND  COORDINATES  OF  OFFSET  NORTH  DIPOLE  AT  ANY 
LATITUDE 

FIRST  GUESS  FIND  OFFSET  NORTH  DIPOLE  AT  DISTANCE  RGC 
ZDEDNP  a  RGC 

100  XODNP  a  XGMED*CPCT  +  YGMED*SP  +  ZD£DNP«CPST 
YODNP  a  -XGMED *SPCT  +  YGMED*CP  -  ZDEDNP*SPST 
ZODNP  a  -XGMED*SP  +  ZDEDNP*CT 

DODNP  a  SQRT (XODNP «XODNP  +  YODNP *YODNP  +  ZODNP*ZODNP ) 

DIFLA  a  DODNP  -  RGC 

IF(  ABS  (DIF LA )  -  1.0E-5)  120,  120,  110 


o  o  o  o  o 


110  ZDEDNP  a  ZDEDNP  -  DIFLA 

4001  FORMAT  (1H  ,  5X, 'ODC  0,  0,  »F7.5,  '  *  GC  X,  Y,  Z  0F'3F8.5, 

1  •  DQDNP  s  'F9.5, '  DIF  OF  'F9.6,  »  AT  LOND  LAT’FIO.4,  F8.4) 

GO  TO  100 
120  CONTINUE 

PHINOF  s  ATAN2(Y0DNP,  XODNP)#RAD 
IF(PHINOF.LT.Q.O)  PHINOF  a  PHINOF  +  360.0 
TNOF  *  -ACCS (ZODNP/DODNP)*RAD  +90.0 

IF(NOPT.EQ.  1  )PRINT  4001, ZDEDNP, XODNP,  YODNP,ZODNP,DODNP, DIFLA, 

1  PHINOF,  TNOF 
SGCLATD  «  SIN (GCLATD/RAD) 

CGCLATD  a  COS  (GCLATD/RAD) 

SGCLOND  =  SIN (GC LOND/ RAD ) 

CGCLOND  a  COS (GCLOND/RAD) 

C  GET  GEOCENTRIC  X  Y  Z  COORDINATES 

XGC  a  RGC*CGCLATD*CGCLOND 
YGC  a  RGC*CGCLATD*SGCLOND 
ZGC  a  RGC*SG^LATD 
GCT  a  (90.0  -  GCLATD)/RAD 
SGCT  a  SIN  (GCT) 

CGCT  a  COS  (GCT) 

FIND  X  Y  Z  IN  LOCAL  COORDINATES  OF  XsO,  Y=0,  Z 

THE  LOCAL  COORDINATE  Z  AXIS  PASSES  THRU  P 
THE  LOCAL  COORDINATE  X,Z  PLANE  CONTAINS  P 
GCR07  a  ATAN2 (YGC, XGC ) 

IF(NOPT. EQ. 1 )PRINT  2001,  XGC,  YGC,  ZGC,  GCROT 
2001  FORMAT ( 1H  ,  4F10.4,  3X,  'XGC,  YGC,  ZGC,  GCROT') 

S GCROT  a  SIN  (GCROT) 

C GCROT  a  COS (GCROT) 

XRL  a  XGC*CGCROT*CGCT  +  YGC*SGCROT»CGCT  -  ZGC*SGCT 
YRL  a  -XGC*S GCROT  +  YGC*CGCROT 

ZRL  a  XGC*CGCROT*SGCT  +  YGC*SGCROT*SGCT  +  ZGC*CGCT 
FORMAT (1H  ,  3F10.4,  13X,  'XRL,  YRL,  ZRL') 

IF(NOPT.EQ.  DPRINT  2002,  XRL,  YRL,  ZRL 

DETERMINE  LOCATION  OF  OFFSET  DIPOLE  CENTER  IN  THESE  SAME 

ROTATED  LOCAL  COORDINATES 

XEDRL  a  XE DF GC*C GCROT *C GCT  +  YEDFGC«SGCROT«CGCT  -  ZEDFGC»SGCT 
YEDRL  a  -XEDFGC*SGCROT  +  YEDFGC*CGCROT 

ZEDRL  a  XEDFGC*CGCROT  *SGCT  +  YEDFGC*SGCROT*SGCT  +  ZEDFGC«CGCT 
IF(NOPT.  EC.  DPRINT  3003,  XEDRL,  YEDRL,  ZEDRL 
3003  FORMAT  (1H  ,  3F10.4,  13X,  'XEDRLM  YEDRL,  ZEDRL') 

XEDRTL  a  XEDRL 
YEDRTL  a  YEDRL 
ZEDRTL  a  ZEDRL  -  ZRL 

IF(NOPT.EQ.  DPRINT  2303,  XEDRTL,  YEDRTL,  ZEDRTL 
2303  FORMAT  (1H  ,  3F10.4,  13X,  'XEDRTL,  YEDRTK,  ZEDRTL') 

C  TRANSUTE  TO  LOCAL  COORDINATE  SYSTEM  WITH  ORIGIN  AT  SURFACE 

XRTL  a  XRL 
YRTL  a  YRL 
ZRTL  a -ZRL 
XEDP  a  XRTL  +  XEDRL 
YE  DP  a  YRTL  +  YEDRL 
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ZEDP  s  ZRTL  +  ZEDRL 

RED  a  SORT  (XEDP*XEDP+YEDP*YEDP+Z  EDP*ZEDP ) 

2302  FORMAT  <1H  ,  3F10.4,  13X,  »XRTL,  YRTL,  ZRTL') 

IF(NOPT.EQ.  DPRINT  2302,  XRTL,  YRTL,  ZRTL 
C.  EARTHS  SURFACE  AT  A  SPECIFIED  ALTITUDE 

C  POSITION  OF  OFFSET  NORTH  DIPOLE  IN  LOCAL  COORDINATE  SYSTEM 

XODNPRa  XO  DN  P  *C  GC  ROT  *C  GCT  +  YODNP«SGCROT*CGCT  -  ZODNP«SGCT 

YODNPRs  -XODNP*SGCROT  +  YODNP^CGCROT 

ZODNPR*  XODNP*CGCROT  *SGCT  +  YODNP  "SGCRGT  *SGCT  +  ZODNP»CGCT 

XODNPT  *  XODNPR 

YODNPT  *  YODNPR 

ZODNPT  *  ZODNPR  -  ZRL 

IF(NOPT,EQ. 1 )PRINT  1103,  XODNPT,  YODNPT,  ZODNPT 
1103  FORMAT (1H  , 10X, 'XODNPT  a 'F 10. 4, 2X, 'YODNPT  *'F10.4,2X, 'ZODNPT  *' 

1  F10.4.3X,  'OFFSET  N  DIPOLE  IN  LQCAL  COORDINATES') 

ROTM  *  ATAN2<Y0DNPT,  XODNPT)  +  PI 

FIND  ANGLE  FROM  GEOGRAPHIC  NORTH 
NEGATIVE  -  ROTATION  FROM  GEOGRAPHIC  NP  CLOCKWISE 
POSITIVE  -  ROTATION  FROM  GEOGRAPHIC  NP  CCW 
SROTM  =  SIN (ROTM) 

CROTM  =  COS (ROTM) 

RCTMD  =  ROTM*RAD 

2327  FORMAT  (1H  ,  F15.5,  3X,  '  ROTM  IN  DEGREES  MEASURED  CCW  SO  -X 

1  WILL  POINT  TOWARD  OFFSET  NORTH  DIPOLE  AXIS') 

IF(NOPT.EQ.  DPRINT  2327,  ROTMD 

C  FIND  COMPONENTS  OF  UNIT  VECTOR  AT  ARBITARY  AZIMUTH  AND  ZENITH 

PLAZ  s  -AZ/RAD  +  PI 
TLZE  *  ZE/RAD 
SPLAZ  =  SIN (PLAZ) 

CPLAZ  *  COS (PLAZ) 

STLZE  s  SIN  (TLZE) 

CTLZE  =  COS (TLZE) 

XLD  x  STLZE*CPLAZ 
YLD  r  STLZE «S PLAZ 
ZLD  a  CTLZE 

IF(NOPT.EQ.  DPRINT  2005,  XLD,  YLD,  ZLD,  AZ,  ZE 
2005  FORMAT  (1H  ,  5F10.4,  3X,  'UNIT  VECTOR  COMPOENTS  AT  AZ  A  ZE ' ) 

C 

C  FIND  COMPONENTS  OF  UNIT  VECTOR  IN  DIPOLE  RADIAL  COORDINATES 

C 

C  ROTATE  AROUND  Y  AXIS  SO  -Z  AXIS  PASSES  THROUGH  XED,  0,  ZED 

C  NEW  VECTOR  IS  VA  a  ZRTL  +  ZEDRTK  +  XEDRTL 

C  ANGLE  BETWEEN  VECTOR  FROM  POINT  LOCAL  ORIGIN  TO  GEOCENTER 

C  AND  VECTOR  FROM  POINT  LOCAL  ORIGIN  TO  XED,  0,  ZED 

C  JIM  LANGWORTHY 'S  FIX 

CAaDPEC  (XEDRTL) 

C  CA  a  ZRTL*ZEDRTL/( ABS (ZRTL)*SQRT (ZEDRTL*ZEDRTL  +  XEDRTL "XEDRTL) ) 

A  a  ACOS(CA) 

IF  (XEDRTL.  GT.  0.0)  A  =  -A 
SA  a  SIN (A  ) 


ADEG  *  A*RAD 

IF(NOPT, EQ, 1 )PRINT  1000,  CA,  A,  SA,  ADEG 

XLP  =  XLD*CA  +  ZLD*SA 

YIP  *  YLD 

ZLP  s-XLD *SA  +  ZLD*CA 

IF(N0PT. EQ. 1 )PRINT  5001,  XLP,  YLP,  ZLP 

5001  FORMAT (1H  ,  3F10.4,  13X,'XLP,  YLP,  ZLP  ’) 

C  ROTATE  AROUND  X  PRIME  AXIS  SO  -Z  PASSES  THROUGH  XED,  YED,  ZED 

CB sDPEC(YEDRTL) 

C  CB  s  Z RTL *ZEDRTL/( ABS ( ZRTL ) *S QRT ( ZEDRTL*ZEDRTL  +  YEDRTL*YEDRTL) ) 

B  ■  ACOS(CB) 

IF (YEDRTL, GT. 0. 0 )  B  =  -B 
SB  =  SIN(B) 

BDEG  =  B*RAD 

IF (NOPT, EQ. 1 )PRINT  1000,  CB,  B,  SB,  BDEG 
XI  PP  s  XLP 

YLPP  s  YLP«CB  +  ZLP«SB 
ZLPP  a -YLP  *SB  +  ZLP*CB 
IF(N0PT. EQ. 1 )PRINT  5002,  XLPP,  YLPP,  ZLPP 

5002  FORMAT (1H  ,  3F10.4,  13X,»XLPP,  YLPP,  ZLPP  •) 

C  ROTATE  AROUND  ZPP  AXIS  SO  -X  AXIS  PASSES  THROUGH  NORTH 

C  OFFSET  DIPOLE  AXIS 

ZLDM  s  ZLPP 

XLDM  =  XLPP«CROTM  +  YLPP«SROTM 
XLDM  «  XLPP*CROTM  -YLPP«SR0TM 
YLDM  s-XLPP*SROTM  +  YLPP«CROTM 
YLDM  =  XLPP«SROTM  +  YLPP«CR0TM 

IF(NOPT. EQ. 1 )PRINT  1101,  XLD,  YLD,  ZLD,  XLDM,  YLDM,  ZLDM 

1101  FORMAT  (1H  ,  'UNIT  VECTOR  IN  LOCAL  COORDINATES  3F8.5,5X, 

1  'UNIT  VECTOR  IN  LOCAL  MAGNETIC  COORDINATES’,  3F10.5) 

C  FIND  AZUMITH  ANGLE  OF  UNIT  VECTOR  IN  LOCAL  DIPOLAR  RADIAL  COOR 

IF( (ABS (YLDM) . GT. SMALL) .OR. (ABS(XLDM) „GT. SMALL))  GO  TO  1102 
PAZMaO.O 
GO  TO  1104 

1102  PAZM  s  ATAN2 (YLDM, XLDM) 

1104  AZM  s  (PI  -  PAZM)*RAD 

IF(AZM.GT. 360.0)  AZM  *  AZM  -  360.0 
ZEM  =  A COS (ZLDM) *R AD 
C  FIND  GAMMA  ANGLE 

GAMMA  a  A COS (YLDM) *R AD 

C  TRANSFORM  TO  OFFSET  DIPOLE  COORDINATES 

XED1  sXGC-XEDFGC 
YED1  =YGC-YEDFGC 
ZEDIaZGC-ZEDFGC 

C  FIND  THE  Z  COORDINATE  IN  OFFSET  DIPOLE  COORDINATES 

ZED2sXED1 *CPST-YED1 «SPST+ZED1 «CT 
C  FIND  THE  GEOMAGNETIC  UTITUDE 

EDLATsRAD*(PI02-AC0S(ZED2/RED) ) 

COSLDA=COS  (EDLAT/RAD) 

STORMERs60.*COSLDA**4./(RED*RED*(1 . +SQRT (1 .-C0SLDA**3.*YLDM) )**2 ) 
RETURN 
END 
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SUBROUTINE  ORBIT (N, TIME, ZLON.ZLAT, RADIUS ) ' 

##*#*«»#*#*####*»#**#•##«###*»###*######*#»# 

THIS  SUBROUTINE  ACCEPTS  INPUT  CONCERNING  SATELLITE 
ORBITS  AND  CALCULATES  THEIR  GEOGRAPHICAL  LOCATION. 

NaOj  ACCEPT  DATA  ON  ORBIT  (FAST  MODE). 

Nsls  ACCEPT  DATA  ON  ORBIT  (COMPLETE  MODE). 

N=2j  CALCULATE  ORBIT  AS  A  FUNCTION  OF  TIME. 

CLEARLY  DATA  MUST  BE  INPUT  BEFORE  COMPUTATIONS. 

ON  DATA  INPUT,  TIME  RETURNS  THE  ORBITAL  PERIOD. 

DURING  ORBIT  CALCULATIONS,  TIME  IS  AN  INPUT  VARIABLE. 

DATA  PI ,  ALTA,  ALTP,  RE ,  E,  A1 ,  A2,  A3,  XI , W1 ,  RMAJ,W2,  FACT 
1, THO, PHO, PSI,XIO  /17*-8000./ 

THESE  VARIABLES  HAVE  BEEN  DATAED  IN  HERE  TO  FORCE  THE  COMPUTER 
TO  CREATE  STORAGE  LOCATIONS  FOR*  THEM  BETWEEN  CALLS  OF  THIS 
SUBROUTINE.  THE  DATAED-IN  VALUES  ARE  NEVER  USED  BY  THE  PROGRAM. 

PI*355./113. 

IF  (N.EQ.2)  GO  TO  1000 

ACCEPT  ORBITAL  DATA. 

ALTAsORBITAL  ALTITUDE  AT  APOGEE  (KILOMETERS). 

ALTPsORBITAL  ALTITUDE  AT  PERIGEE  (KILOMETERS). 

RE sRADIUS  OF  EARTH  (KILOMETERS). 

E=ORBIT  ECCENTRICITY. 

A1 sORBITAL  INCLINATION  (DEGREES). 

A2aINITIAL  LONGITUDE  OF  ASCENDING  NODE  (DEGREES). 

A3*INITIAL  DISPLACEMENT  FROM  ASCENDING  NODE  (DEGREES). 
XIsDISPLACEMENT  OF  PERIGEE  FROM  ASCENDING  NODE  (DEGREES). 

WHAT  IS  THE  ALTITUEDE  AT  APOGEE? 

WRITE (6, 420) 

ACCEPT  »,  ALTA 

WHAT  IS  THE  ALTITUDE  AT  PERIGEE? 

WRITE (6, 400) 

ACCEPT  »,  ALTP 
RE=6371.2 

E  s (ALTA -ALTP ) / ( ALTA+ALTP+2 . *RE ) 

IF  (E.LT.. 00001 )  EsO. 

WHAT  IS  THE  ORBITAL  INCLINATION? 

WRITE  (6, 405) 

ACCEPT  «,  A1 
A2s0. 

A3=0. 


OOO  OOO  O  O  O  O  O  O  O  O  o  O  o  o  o 


c 

c 


100 


400 

405 

410 

1 

415 

1 

420 

425 


1000 


XI  =0. 

IF  (N.EQ.0)  GO  TO  100 

WHAT  IS  THE  INITIAL  LONGITUDE  OF  THE  ASCENDING  NODE? 

WRITE  (6 ,410) 

ACCEPT  »,  A2 

WHAT  IS  THE  INITIAL  DISPLACEMENT  FROM  THE  ASCENDING  NODE? 

WRITE(6, 415) 

ACCEPT  »,  A3 
IF  (E.EQ.O.)  GOTO  100 

WHAT  IS  THE  DISPLACEMENT  OF  THE  PERIGEE  FROM  THE  ASCENDING  NODE? 

WRITE (6, 425) 

ACCEPT  »,  XI 
CONTINUE 

W1 sANGULAR  VELOCITY  OF  EARTH  (RADIANS/SEC). 

RMAJsSEMI -MAJOR  AXIS  (KILOMETERS). 

W2=MEAN  ORBITAL  ANGULAR  VELOCITY  (RADIANS/ SECOND) . 

TIME=ORBITAL  PERIOD  (SECONDS). 

FACTsA  USEFUL  FACTOR. 


W1-7. 27E-5 

RMAJ=(ALTP+RE)/(  1.  -E ) 

W2=1 . 24E-3*(RE/RMAJ  )»«1 . 5 

TIME=2.«PI/W2 

FACT=SQRT((1.4E)/(1.-E)) 

DEFINE  MORE  USEFUL  ANGLES. 

THO=PI*A1/1  80. 

PHO*PI»(A2-90.)/180. 

PSIsPI*A3/1  80. 

XIO=PI«n/180. 

FORMAT  STATEMENTS. 

FORMAT (IX, ’ENTER  ALTITUDE  AT  PERIGEE  (KILOMETERS):  ’$) 
FORMAT (IX, 'ENTER  ORBITAL  INCLINATION  (DEGREES):  »$ ) 
FORMAT (IX, 'ENTER  INITIAL  LONGITUDE  OF  ASCENDING  NODE', 

’  (DEGREES):  '$) 

FORMAT (IX, 'ENTER  INITIAL  DISPLACEMENT  FROM  ASCENDING', 

'  NODE  (DEGREES  ):  '$) 

FORMAT (IX, 'ENTER  ALTITUDE  AT  APOGEE  (KILOMETERS):  '$  ) 
FORMAT (IX, 'ENTER  DISPLACEMENT  OF  PERIGEE  FROM', 

»  ASCENDING  NODE  (DEGREES):  '$ ) 

RETURN 

CONTINUE 
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COMPUTE  SATELLITE  POSITION. 

QM=MEAN  ANOMALY. 

IF  (E.NE.O.)  GOTO  1009 
QM=W2«TIME+PSI 
GOTO  1010 

1009  CONTINUE 

YSS  =  (PS I -XI 0 )  /2 . 

QE0 s2 .* ATAN2 (SIN ( YSS ) , FACT»C0S ( YSS ) ) 
QM0sQE0-E*SIN(QE0) 

CMsW2»TIME-QM0 

1010  CONTINUE 

QE*EC CENTRIC  ANOMALY. 

QE=QM 

DEL=1. 

1014  CONTINUE 

QTEMP  =QE 
QE=QM+E«SIN(QE) 

DEL =QE -QTEMP 

IF  ( ABS (DEL ) . GT . „ 000 1 )  GOTO  1014 

QT=TRUE  ANOMALY. 

IF  (E.NE.O.)  GOTO  1019 
QT=QE 
GOTO  1020 

1019  CONTINUE 

QECYC=INT(QE/2./PI ) 

QE  RED  sQE-2 .  *PI  *QEC  YC 
AA 1  sFACT*SIN (QERED/2 . ) 

AA2  sCOS  (QERED/2.) 

QTREDr2.*ATAN2(AA1,AA2) 

IF  (QTRED.LT.O.)  QTRED  =QTRED+2 . *PI 
QT=QTRED+2.*PI  *QECYC+XIO 

NOTE:  ANOMALY  COMPUTATIONS  ARE  DONE  FROM  PERIGEE 
WHILE  ORBIT  COMPUTATIONS  BELOW  ARE  DONE  FROM  THE 
ASCENDING  NODE.  THE  FACTOR  OF  XIO  CORRECTS  THIS. 

1020  CONTINUE 

ZLATsLATITUDE. 

R1«SIN  (THO  )*SIN(QT ) 

THP  =ACOS (R 1 ) 

ZLAT*90.-1  80.*THP/PI 

ZLONsLONGITUDE. 


Cl  o  o 


C 


RP=, 5*(PI/2, +THO ) 

RM=.5*  (PI/2, -THO ) 

RFs,  5*(PI/2,~QT) 

IF  (SIN (RF).NE. 0, )  GOTO  1029 
PHP=PI+PH0-W1  *TIME 
GOTO  1030 

1029  CONTINUE 

S1*SIN(RM)*COS(RF)/SIN(RP)/SIN(RF> 
S2 sCOS (RM)»COS (RF)/COS (RP)/SIN (RF) 
SUM  sATAN  (S 1  )+ATAN  (S2) 

PHP=PHO-W1  *TIMF+SUM 

1030  CONTINUE 

IF  (PHP.GE.O. )  GOTO  1034 
PHP=PHP+2.«PI 
GOTO  1030 

1034  CONTINUE 

IF  (PHP.LT.2. *PI )  GOTO  1035 
PHP=PHP-2.*PI 
GOTO  1034 

1035  CONTINUE 
Z LON  si  80.»PHP/PI 

RADIUS sALTITUDE  (KILOMETERS). 

RADIUS  =RMAJ*( 1 . -E  *COS (0 i  )  )-R£ 

RETURN 
END 
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PROGRAM  LET 

THIS  PROGRAM  COMPUTES  INTEGRAL  AND  DIFFERENTIAL  LET  SPECTRA. 

THE  PARTICLE  SPECTRA  IN  THE  INTER  PLANETARY  MEDIUM  ARE  IN 
THE  CRF  SUBROUTINE. 

THE  GEOMAGNETIC  CUTOFF  IS  APPLIED  BY  THE  CUT  SUBROUTINE. 
PROPAGATION  THROUGH  SHIELDING  IS  HANDLED  IN  THE  INSIDE  ROUTINE. 
THE  DEDXSI  ROUTINE  RETURNS  LET  IN  SILICON. 

THE  CONTROL  PARAMETERS  ARE: 

Y  a  YEAR  OF  THE  PARTICLE  ENVIRONMENT  O  975.  144  WAS  SOLAR  MIN. 
AND  1980.598  WAS  SOLAR  MAX.); 

IT  a  ORBITAL  INDEX  (0  a  HIGH  ORBIT:  NO  OUXOf F,  NO  TRAPPED 
PROTONS;  1  a  LCW  ORBIT:  CUTOFF  AND  TRAPPED  PROTONS); 

MQ  a  INTERPLANETARY  MEDIUM  WEATHER  INDEX  (SEE  NOTES  IN  CRF 
ROUTINE); 

THK  =  SHIELDING  THICKNESS  .1  ALUMINUM  EQUIVALENT)  IN  G/CM**2; 


C  IDEM  a  NO.  OF  LOG  STEFS  IN  ENERGY; 

C  IELM  a  LIGHTEST  ELEMENT  INCLUDED; 

C  JELM  a  HEAVIEST  ELEMENT  INCLUDED  (JELM  MUST  NOT  EXCEED  92). 

C  THE  OUTPUTS  ARE  THE  DIFFERENTIAL  AND  INTEGRAL  LET  SPECTRA. 

C  THE  UNITS  ARE  PARTICLES/((M**2)*STER*SEC*MEV*CM**2/G  )  AND 

C  PARTICLES/((M**2)*STER*SEC)  RESPECTIVELY 

C  VERSUS  MEV»CM*«2/G. 


C 


C  DXMAX  AND  DXMIN  THE  MAXIMUM  AND  MINIMUM  LET  VALUES  CONSIDERED, 

r  LMIH  IS  THE  MINIMUM  ENERGY  CONSIDERED. 


C  ABOVE  EMAX,  LET  IS  CONSTANT. 

C 

REAL  NO 01), INSIDE 

CHARACTER *12  ALPHA, BETA, GAMMA 

DIMENSION  DX(101  ),F(  1000), DEDX(1 000) 

DATA  DXMAX /I, 1E+5/, DXMIN /I. 6/, EMIN/. 1 / , EMAX/20000,/ 
DATA  F/1 000*0.0 /.GAMMA/’  ' /,NDIF/0/,NINT/0/ 


C  ASK  FOR  ALL  THE  CONTROL  DATA  NEEDED  FOR  A  RUN. 

C 

WRITE  (6,101) 

101  FORMAT ( '  ENTER  THE  NAME  OF  THE  DIFFERENTIAL  LET  SPECTRUM'/ 

1  '  (ENTER  NOTHING  IF  YOU  DO  NOT  WANT  ONE):  ',$) 

ACCEPT  102,  ALPHA 


102 


103 


104 

105 


FORMAT (A 12) 

IF( ALPHA. EQ. GAMMA)  GO  TO  104 
OPEN (UNIT a«0, FILE aALPHA) 

WRITE  (6,103) 

FORMAT ('  ENTER  THE  NUMBER  OF  POINTS  IN  THE  TABULATION  OF'/ 
'  THE  DIFFERENTIAL  LET  SPECTRUM  (MAXIMUM:  1000):  ',$) 

ACCEPT  *,NDIF 
WRITE  (6,105) 

FORMAT (’  ENTER  THE  NAME  OF  THE  INTEGRAL  LET  SPECTRUM’/ 

'  (ENTER  NOTHING  IF  YOU  DO  NOT  WANT  ONE):  »,$) 

ACCEPT  102,  BETA 

IF( BETA. EQ. GAMMA)  GO  TO  107 

OPEN  (UNITa45 ,  FILE  aBETA ) 
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WRITE  (6,106) 

106  FORMAT ( '  ENTER  THE  NUMBER  OF  POINTS  IN  THE  TABULATION  OF'/ 

1  '  THE  INTEGRAL  LET  SPECTRUM  (MAXIMUM:  1000):  '  ,$) 

ACCEPT  #,NINT 

107  NPT  =MAXO (NDIF, NINT ) 

WRITE  (6,108) 

108  FORMAT ('  ENTER  THE  ATOMIC  NUMBER  OF  THE  LIGHTEST  ELEMENT'/ 

1  '  TO  BE  INCLUDED  IN  THE  LET  SPECTRA:  ',$) 

ACCEPT  *,IELM 
WRITE  (6,109) 

109  FORMAT ('  ENTER  THE  ATOMIC  NUMBER  (<93)  OF  THE  HEAVIEST  ELEMENT'/  * 

1  '  TO  BE  INCLUDED  IN  THE  LET  SPECTRA:  ',$) 

ACCEPT  *, JELM 
IF  (JELM. GT. 92)  JELM=92 
WRITE  (6,110) 

110  FORMAT ('  ENTER  THE  YEAR  FOR  WHICH  YOU  WANT  THE  SPECTRUM'/ 

1  ’  (1975.144  FOR  SOLAR  MIN.}  1980.598  FOR  SOLAR  MAX.):  ',$) 

ACCEPT  *,Y 
WRITE  (6,111) 

111  FORMAT ('  ENTER  THE  ORBITAL  INDEX:  0  FOR  NO  CUTOFF  AND  NO  ', 

1  TRAPPED  PROTONS}'/*  1  FOR  CUTOFF  AND  TRAPPED  PROTONS', 

2'  (NOTE:  STASS.D/tT  AND  .7 TRANS.  DAT '/ 

3’  MUST  BE  SUPPLIED  WHEN  1  IS  ENTERED  HERE):  ',$) 

ACC,;  yf  5, it 
WRITE  (6,112) 

112  FORMAT ('  ENTEi-:  THE  INVERPUNETJU \  WEATHER  INDEX:  ',$) 

ACCEPT  *,  MQ 

WRITE  (6,113) 

113  FORMAT ('  ENTER  THE  , SHIS., DING  THICKNESS  IN  INCHES  OF'/ 

1  '  ALUMINUM  (OR  EQUIVA  LENT ):  ',$) 

ACCEPT  *,THK 

CONVERT  V.  S  THICK  ESS  TO  G/CM**2  (ALUMINUM  EQUIVALENT). 

THK*THK«6.  i.528 
WRITE  (6,11:')  IELM, JELM 

116  FORMAT(//'  THE  LET  SPECTRA  WILL  INCLUDE  THE  CONTRIBUTIONS'/ 

1  '  OF  THE  ELEMENTS  ',12,'  THROUGH  ',12,'.') 

IF(NDIF.GT.O)  WRITE  (6,114)  ALPHA, NPT 
IF(NINT.GT.O)  WRITE  (6,115)  BETA, NPT 

114  FORMAT ( '  THE  FILE  \A12,'  WILL  CONTAIN  THE  DIFFERENTIAL  LET'/ 

1  '  SPECTRUM,  TABULATED  IN  ',14,'  DATA  POINTS.') 

115  FORMAT ( '  THE  FILE  ' ,A12, *  WILL  CONTAIN  THE  INTEGRAL  LET'/ 

1  '  SPECTRUM,  TABULATED  IN  ',14,'  DATA  POINTS.') 

WRITE  (6,117)  Y 

117  FORMAT ( '  THESE  LET  SPECTRA  WILL  BE  FOR  THE  YEAR  ' ,F8. 3. ' > 

IF(IT.EQ.O)  WRITE  (6,118) 

118  FORMAT ('  THESE  LET  SPECTRA  WILL  BE  UNSHIELDED  BY  THE'/ 

1  '  EARTH  "S  MAGNETIC  FIELD,  AND  THE  CONTRIBUTION  FROM  TRAPPED'/ 

2  '  PROTONS  WILL  NOT  BE  INCLUDED. ' ) 

IF(IT.EQ.I)  WRITE  (6,119) 

119  FORMAT ('  THESE  LET  SPECTRA  WILL  INCLUDE  THE  EFFECTS  OF'/ 
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1  •  GEOMAGNETIC  SHIELDING  ACCORDING  TO  THE  GEOMAGNETIC  CUTOFF'/ 

2  '  TRANSMITTANCE  FUNCTION  TABULATED  IN  GTRANS.DAT.  THE'/ 

3  '  CONTRIBUTION  FROM  TRAPPED  PROTONS  WILL  BE  INCLUDED'/ 

4  '  ACCORDING  TO  THE  ORBIT -AVERAGED  TRAPPED  PROTON  SPECTRA*/ 

5  '  TABULATED  IN  STASS.DAT. ' ) 

WRITE  (6,120)  MQ 

120  FORMAT ( '  THE  INTERPLANETARY  WEATHER  INDEX  IS  ',12,'.') 

WRITE  (6,121  )  THK 

121  FORMAT ( '  THE  SHIELDING  THICKNESS  IS  \F7.3,'  GM/CM»*2 
1  ALUMINUM  (OR  EQUIVALENT).'//) 

4  DLLsALOG  (DXMIN) 

DD=(ALOG  (DXMAX  )-DLL)  /NPT 
DLXsDLL 
IDEMsIOO 
JDEMsIDEM+1 

SET  UP  LOG  ENERGY  STEPS. 

ELNXsALOG(EMIN) 

DELNa ( ALOG (EMAX )-ELNX ) /IDEM 

DO  ONE  ELEMENT  AT  A  TIME. 


DO  41  IZ=IELM, JELM 

SET  LOWER  BIN  BOUNDARIES. 

L»1 

E sEMIN 
Ksl 

ELN=ELNX 

ZsFLOAT(IZ) 

GET  LET. 

DX ( 1 )sDEDXSI (Z, E) 

GET  INTERNAL  FLUX. 

FLsINSIDE (Z,E, Y,MQ, IT, THK) 

RETURN  TO  HERE  AFTER  FLUX  IN  EACH  BIN  HAS  BEEN  CALCULATED. 

10  LsL+1 

RETURN  TO  HERE  IF  D(LET  )/DE  IS  SINGULAR  TO  INCREASE  THE  BIN 
WIDTH. 

11  KsK+1 


WHEN  K  EXCEEDS  JDEM  YOU  ARE  THROUGH. 
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IF(K.GT.JDEM)  GO  TO  20 

STEP  TO  UPPER  ENERGY  BIN  BOUNDARY. 

ELNsELN+DELN 

COMPUTE  UPPER  ENERGY  BIN  BOUNDARY. 

EUsEXP(ELN) 

COMPUTE  UPPER  LET  BIN  BOUNDARY.  * 

DUsDEDXSI (Z,EU) 

TEST  FOR  D(LET )/DE  SINGULARITY. 

T=ABS(DX(L-1  )-DU)/DU 
IF(T.LT.  1.0E-6)  GO  TO  11 

COMPUTE  FLUX  IN  THE  BIN. 

FUsINSIDE (Z,EU, Y,MQ, IT,THK) 

N  (L-1  )s0.5*(FU4FL  )*  (EU-E  ) 

THE  UPPER  BIN  BOUNDARY  BECOMES  THE  LOWER  BOUNDARY  OF  NEXT  BIN. 

DX(L)*DU 

EsEU 

FL=FU 

CLOSE  THE  LOOP. 

GO  TO  10 

DO  THE  HIGH  ENERGY  TAIL  IN  THE  LAST  BIN,  ASSUMING  AN  E»«(-2.5) 

POWER  LAW  FOR  HIGH  ENERGY  COSMIC  RAYS. 

LL=L-2 

N (LL+1 )s0.6666#FL*E 

PRESERVE  THE  LET  BIN  BOUNDARIES. 

DLU=DLL 

DLsDXMIN 

RE-BIN  EACH  ELEMENT  IN  A  COMMON  SET  OF  LET  BINS. 

DO  41  KKsI ,  NPT 

SET  UP  THE  COMMON  LET  BIN  BOUNDARIES. 

DLUsDLU+DD 
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DUsEXP(DLU) 
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C  SPIN  THROUGH  THE  ELEMENT  BIN  SET  ADDING  FLUX  TO  THE 

C  COMMON  BIN. 

C 

DO  40  Isl.LL 
DX1  =DX  (I ) 

DX2=DX(I+1 ) 

C 

C  SKIP  THE  LOOP  FOR  ELEMENT  BINS  NOT  OVERLAPPING  THE  COMMON  BIN. 

C 

IF  ( (DX1 ,GE. DU ) .AND. (DX2. GE. DU ) )  GO  TO  40 
IF  ( (DX1.LE.DL) .AND. (DX2. LE.DL ) )  GO  TO  40 
IFCDX1.LT. DX2)  GO  TO  35 
IFCDX1.GE.DU)  GO  TO  32 
IFCDX2.LT.DL)  GO  TO  31 
FCKK)sF(KK)+N  (I ) 

GO  TO  40 

31  FCKK  )sp(KK  )+N  (I  )»CDX1-DL)/(DX1-DX2) 

GO  TO  40 

32  IFCDX2.GE.DL)  GO  TO  33 
FCKK)=FCKK)+N  (I )» CDU-DL)  /  (DX1 -DX2  ) 

GO  TO  40 

33  FCKK  )=F(KK  )+N  Cl  )»CDU-DX2  )/CDX1-DX2  ) 

GO  TO  40 

35  IFCDX2.GE.DU)  GO  TO  37 
IFCDX1.LT.DL)  GO  TO  36 
FCKK  )=F(KK  )+N (I ) 

GO  TO  40 

36  FCKK  )*FCKK  )+N  (I  )*CDX2-0L  )/CDX2-DX1  ) 

GO  TO  40 

37  IFCDX1.GE.DL)  GO  TO  38 

FCKK  )=FCKK)+N  (I )« CDU-DL) /CDX2-DX1  ) 

GO  TO  40 

38  FCKK)=FCKK)+N(I  )«CDU-DX1  )/CDX2-DX1  ) 

40  CONTINUE 
C 

C  DO  THE  INTEGRAL  BIN  FOR  THE  HIGH  ENERGY  TAIL. 

C 

IF( (DX2.GE. DL)  .AND. (DX2.  LE. DU ) )  FCKK  )*F(KK  )+N  (LL+1  ) 

C 

C  FINISH  STEPPING  THROUGH  THE  COMMON  LET  BINS 

C  THEN  GO  BACK  AND  DO  THE  NEXT  ELEMENT. 

C 

41  DL=DU 
C 

C  COMPUTE  THE  DIFFERENTIAL  LET  SPECTRA. 

C 

DL=DXMIN 
DO  50  K= 1 , NPT 
DLX=DLX+DD 
DU^EXP  (DLX ) 


87 


G  a  F  (K)/(DU-DL  ) 

DEDX(K)  a  0.5*<DU+DL) 

IF(NDIF.GT.O)  WRITE  (40,600  )  G,  DEDX(K) 
FORMAT (2X, E12. 5, 2Xf  E12. 5 ) 

Dt  =DU 
SUMaO.O 

INTEGRATE  THE  LET  SPECTRA. 

DO  60  IaNPT,1,-1 
SUMsSUM+F(I) 

DLaDEDX(I) 

IF(NINT.GT.O)  WRITE (45, 600 )SUM,DL 

CONTINUE 

END 


FUNCTION  DED)CSI(Z,E) 

C  THIS  SUBROUTINE  RETURNS  THE  STOPPING  POWER  IN  SILICON  IN 

C  MEV«CM»*2/G. 

C  Z  IS  THE  ION  ATOMIC  NUMBER. 

C  E  IS  THE  ION  ENERGY  IN  MEV/AMU. 

C  ZN  ARE  THE  ATOMIC  NUMBERS  FOR  WHICH  THE  STOPPING  POWERS  HAVE 

C  BEEN  TABULATED. 

C  KK  ASSOCIATES  EACH  ELEMENT  WITH  ONE  WHOSE  STOPPING  POWER  HAS 

C  BEEN  TABULATED. 

C 

DIMENSION  EN  (77 ) , RN (9 , 77 ) , ZN (9 ) , KK (92 ) 

DATA  IST/0/,ZN/1.,2.,6.,8.,18.,26.,36.,54.,S2/ 

DATA  IZS, J/0,3/,KK/1,2*2,3«3,7*4,8*5,9»6 
1  ,14*7,28*8, 20*9/ 

C 

C  THE  STOPPING  POWER  TABLES  ARE  IN  THE  SILICN.DT1  FILE. 

C 

C 

C  READ  THE  STOPPING  POWER  TABLES  UPON  THE  FIRST  ENTRY. 

C 

IF(IST.NE.O)  GO  TO  5 

1ST  -i 

OPEN (UNIT=5 1 , READONLY, SHARED, STATUS  a ’OLD » ,FILE='SILICN. DTI') 

DO  3  1*1,77 

READ (5 1,2)  EN(I),(RN(J,I),Ja1,9) 

2  FORMAT  (6X,F9.  3,9  (3X,E9.3)> 

3  CONTINUE 
CLOSE  (UNIT  =51) 

5  CONTINUE 

IZ=Z+0.5 
C 

C  COMPUTE  THE  APPROXIMATE  STARTING  INDEX  FOR  TABLE  LOOKUP. 

C 

I J= (ALOG (E )+4 . 065 )*4 . 85567-1 . 0 
IJ  =MAXO (I J, 1  ) 

IJ=MIN0(IJ,77> 

C 

C  SEARCH  FOR  THE  ENERGY  TABULAR  POINT  JUST  ABOVE  E. 

C 

10  CONTINUE 

DO  11  I  si  J,  77 
IF(EN(I).LT.E)  GO  TO  11 
J=I 

GO  TO  13 

1 1  CONTINUE 
J=77 

C 

C  CALCULATE  THE  INTERPOLATION  FACTOR. 

C 

13  CONTINUE 

F=(E-EN(J-1  ))/(EN(J  )-EN  (J-1  ) ) 
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FsAMINKF,  1.0) 

FIND  THE  NEAREST  TABULATED  ELEMENT. 

14  CONTINUE 
KsKK(IZ) 

INTERPOLATE  THE  STOPPING  POWER  FOR  THE  TABULATED  ELEMENT. 

DEDXSIsF*  (RN  (K,  J  )-RN  (K ,  J-1  )  )+RN  (K ,  J-1 ) 

ZNEsZN(K) 


IF  THE  ION  IS  THE  ONE  TABULATED,  YOU  ARE  DONE. 

IF  (Z.EQ.ZNE )RETURN 
ZE=Z 

ECa1.+2.«Z 

IF  THE  ENERGY  IS  HIGH  ENOUGH,  THE  ION  IS  FULLY  STRIPPED. 

IF(E.GT.EC)  GO  TO  21 

COMPUTE  THE  EQUILIBRIUM  CHARGE  STATE. 

BaSQRTd  .-(E«. 001073927+1  )**(-2. ) ) 
ZEsZ*(1.-EXP(-130.»B»(Z**(-.6666)))) 

ZNEsZNE*(1 .  -EXP  C-130.*B*(ZNE**(-.  6666) ) ) ) 

INTERPOLATE  TO  THE  ION  IN  QUESTION. 

21  CONTINUE 

DEDXSI=(ZE«ZE)/(ZNE«ZNE)*DEDXSI 

RETURN 

END 


FUNCTION  INSIDE  < Z, E, Y,M, IT,TH ) 


t 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 
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TRANSPORTS  THE  PARTICLE  SPECTRA  THROUGH  THE  SPACECRAFT  WALL 
CONSIDERING  ENERGY  LOSS  EFFECTS  AND  NUCLEAR  FRAGMENTATION, 

BUT  DOES  NOT  KEEP  TRACK  OF  THE  FRAGMENTS. 

Z  IS  THE  ATOMIC  NUMBER  OF  THE  ION. 

E  IS  THE  ENERGY  IN  MEV/AMU. 

Y  IS  THE  YEAR  OF  THE  ENVIRONMENT. 

M  IS  THE  INTERPLANETARY  MEDIUM  WEATHER  INDEX  (SEE  CRF.FOR ) . 

IT  IS  THE  ORBITAL  INDEX:  1  FOR  LCW  ORBITS,  0  FOR  HIGH  ORBITS. 
TH  IS  THE  THICKNESS  OF  ALUMINUM  SHIELDING  IN  G/CM**2. 

THIS  ROUTINE  RETURNS  THE  INTERNAL  FLUX  IN  PARTICLES/( (M««2 )* 
STER»SEC*MEV/U). 

AN  ARE  THE  AVERAGE  ION  MASSES  EOR  EACH  ELEMENT  (TAKEN  FROM 
THE  HANDBOOK  OF  CHEMISTRY  AND  PHYSICS,  59TH  EDITION). 

REAL  INSIDE 
COMMON/MASS/ AN (92 ) 

DATA  (AN(I),I=1,92)/1.,4.  ,6.9, 

1  9. ,  1 0. 8, 12. ,  1 4. ,  1 6. ,  1  9. , 20. 2,23. , 24.  3,27.  ,28. , 

2  31., 32., 35. 5, 39. 9, 39., 40., 45., 47. 9, 50. 9, 

3  52., 54. 9, 55.  8, 58.  9, 58. 7, 63. 5, 65.  4, 69. 7,72. 6, 

4  74.9, 79.  ,79. 9, 83. 8, 85.5, 87. 6,  88. 9,  91 . 2, 92. 9, 95. 9, 97. ,  101 . , 

5  102.9,106.4,107.9,112.4,114.8,118.7,121.8,127.6,126.9,131.3, 

6  132.9, 137.3, 138.9, 140. 1, 140.9, 144.2, 145. ,150.4, 152. ,157.3, 

7  158. 9, 162. 5, 164. 9, 167. 3, 168. 9, 173.,  175.,  178. 5, 180. 9, 183. 9, 

8  186.2,  190.2,  192.2,  195.1,  197. ,200.6,204.4,207.2,209. ,209. , 

9  210., 222., 223., 226., 227., 232., 231., 238./ 

IZ=Z+.  5 
AsAN(IZ) 

COMPUTE  THE  TOTAL  INELASTIC  CROSS  SECTION  USING  C.H. 

TSAO'S  VERSION  OF  THE  BRANDT -PETERS  FORMULA. 

USE  A  SIMPLER  FORMULA  FOR  PROTONS. 

IF  (Z.NE.1.)  COEF  =  (  (6, Q2E+23  )*(5.  OE-26) 

1  «(A»»(1./3.)+AN(13)**(1./3.)-0.4)«»2)/AN(15) 

IF  (  Z.EQ.  1. )  COEFsO. 

FIND  THE  RESIDUAL  RANGE  OF  THE  ION  INSIDE  THE  SPACECRAFT. 
RsRAL(Z,A,E) 

FIND  THE  RESIDUAL  RANGE  OF  THE  SAME  ION  OUTSIDE  THE  SPACECRAFT, 


R=R+TH 


FIND  THE  CORRESPONDING  ENERGY  OUTSIDE. 
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ES=EAL(Z,  A,R ) 

FIND  THE  STOPPING  POWER  INSIDE  AND  OUTSIDE. 

SS =DEDXAL(Z»  ES ) 

S=DEDXAL( Z,E ) 

GET  THE  FLUX  OUTSIDE  AND  COMPUTE  IT  INSIDE. 

ITsO  FOR  AN  ORBIT  SO  HIGH  THAT  CUTOFF  IS  INSIGNIFICANT 
AND  TRAPPED  PROTONS  ARE  INSIGNFICANT. 

ITsI  FOR  LOWER  ORBITS:  INCLUDE  CUTOFF  AND  TRAPPED  PROTONS. 

IF  (IT.EQ.O)  INSIDE=CRF(Z,ES,Y?M)«(SS/S)»EXP(-COEF»TH) 

IF  (IT.GT.O  )  INSIDE  =CUT(Z,ES,YtM)*(SS/S)*EXP  (-COEF*TH  ) 

RETURN 

END 


FUNCTION  DEDXAL(Z,E) 

C 

C  THIS  SUBROUTINE  RETURNS  THE  STOPPING  POWER  IN  ALUMINUM  IN 

C  MEV»CM**2/G. 

C  Z  IS  THE  ION  ATOMIC  NUMBER. 

C  E  IS  THE  ION  ENERGY  IN  MEV/AMU. 

C 

C  ZN  ARE  THE  ATOMIC  NUMBERS  FOR  WHICH  THE  STOPPING  POWER  HAS  BEEN 

C  TABULATED. 

C  KK  ASSOCIATES  EACH  ELEMENT  WITH  ONE  WHOSE  STOPPING  POWER 

C  HAS  BEEN  TABULATED. 

C 

DIMENSION  EN (66 ) , RN (9 , 66 ) , ZN (9  > , KK (92 ) 

DATA  IST/0/,ZN/1.,2.,6.,8.,18.,26.,36.  ,54.,92/ 

DATA  IZS,  J/0 , 3/  ,KK/1 , 2*2 , 3*3 , 7*4.  8*5 , 9*6 
1  ,14*7, 28*8, 20*9/ 

C 

C  THE  STOPPING  POWER  TABLES  ARE  ON  THE  ALUM NM. DTI  FILE. 

C 

C 

C  READ  THE  STOPPING  POWER  TABLES  UPON  THE  FIRST  ENTRY. 

C 

IF(IST.NE.O)  GO  TO  5 
1ST  si 

OPEN  (UNITsSO ,  READONLY,  SHARED,  STATUS s  'OLD ' , FILEs  ’ALUMNM .  DT 1 » ) 

DO  3  1=1,66 

READ (50, 2)  EN(I),(RN(J,I),J=1,9) 

2  FORMAT (6X,F9. 3,9 (3X.E9.3)) 

3  CONTINUE 

5  CONTINUE 

IZ=Z+0. 5 
C 

C  COMPUTE  THE  APPROXIMATE  STARTING  INDEX  FOR  TABLE  LOOKUP. 

C 

I J  s ( ALOG (E )+4 . 605 ) *4 . 85567-1 . 0 
IJ  sMAXO (I J, 1  ) 

I  J=MINO(I  J,  66) 

C 

C  SEARCH  FOR  THE  ENERGY  TABULAR  POINT  JUST  ABOVE  E. 

C 

10  CONTINUE 

DO  11  I  si J, 66 
IF(EN(I).LT.E)  GO  TO  11 
J=I 

GO  TO  13 

1 1  CONTINUE 
Js66 

C 

C  CALCULATE  THE  INTERPOLATION  FACTOR. 

C 

13  CONTINUE 
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F=  (E-EM  (J-1 )  )/(EN  (J  )-EN  (J-1  ) ) 

FaAMINKF,  1.0) 

FIND  THE  NEAREST  TABULATED  ELEMENT. 

14  CONTINUE 

KsKK(IZ) 

INTERPOLATE  THE  STOPPING  POWER  FOR  THE  TABULATED  ELEMENT. 

20  CONTINUE 

DEDXALaF*  (RN(K,J)-RN(K,J-1  >)+RN(K,  J-1  ) 

ZNEsZN(K) 

IF  THE  ION  IS  THE  ONE  TABULATED^  YOU  ARE  DONE. 

IF  (Z.EQ.ZNE)  RETURN 
ZE=Z 

ECsl  . +2.*Z 

IF  THE  ENERGY  IS  HIGH  ENOUGH,  THE  ION  IS  FULLY  STRIPPED. 

IF(E.GT.EC)  GO  TO  21 

COMPUTE  THE  EQUILIBRUM  CHARGE  STATE. 

B=SQRT(1  .-(E«. 001073927+1  )*»(-2. ) ) 

ZE=Z*(  1.  -EXP (-130.*B*( Z**(-. 6666) ) ) ) 

ZNE=ZNE«(1.-EXP  (-130.  *B»(ZNE««<-.  6666)))) 

INTERPOLATE  TO  THE  ION  IN  QUESTION. 

21  CONTINUE 

DEDXALs (ZE*ZE )/( ZNE*ZNE ) *DEDXAL 

RETURN 

END 
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FUNCTION  RAL(Z,A,E) 

C  THIS  ROUTINE  RETURNS  THE  RESIDUAL  RANGE  IN  ALUMINUM  IN  G/CM**2. 

C  Z  IS  THE  ION  ATOMIC  NUMBER. 

C  E  IS  THE  ION  ENERGY  IN  MEV/AMU. 

C  A  IS  THE  ION  ATOMIC  MASS. 

ZN  ARE  THE  ATOMIC  NUMBERS  OF  THE  TABULATED  ELEMENTS. 

KK  ASSOCIATES  EACH  ELEMENT  WITH  A  TABULATED  ELEMENT. 

AN  ARE  THE  MEAN  ATOMIC  MASSES  OF  THE  TABULATED  ELEMENTS. 

DEDX  ARE  STOPPING  POWER  DATA  IN  (MEV/AMU )/ (G/CM**2 ) . 

DIMENSION  EN (66 ) ,  RN (9 , 66 ) ,  ZN (9  ) , AN (9  )  ,KK (92 ) , DEDX (9 ) 

DATA  IST/O/  ,ZN/1. ,  2. ,  6.  f  8. ,  1  8. ,  26. ,  36.  ,54. ,  92  / 

DATA  IZS,  J/0 , 3/.KK/1 , 2*2 , 3*3 , 7*4 , 6*5 , 9*6 
1  ,14*7,28*8,20*9/ 

DATA  AN/1 . 008, 4. 0026, 12.01115,15. 99994, 39. 948, 55. 847 
1  ,83.8,131.3,238.03/ 

DATA  DEDX/1. 706, 1.720, 5. 174, 6. 913, 14. 12, 21. 22, 27. 34 
1  ,40.05,68.69/ 

SET  ENTRY  FLAG  FOR  RANGE. 

NTRY=0 

READ  IN  THE  TABLES  UPON  THE  FIRST  ENTRY. 

IF(IST.NE.O)  GO  TO  5 

FILE  ALUMNM.DT2  CONTAINS  THE  RANGE-ENERGY  TABLES. 

1  CONTINUE 

OPEN (UNIT=60, READONLY, SHARED, STATUS* ’OLD • ,FILEs’ALUMNM. DT2  • ) 
1ST  si 

DO  3  1*1,66 

READ(60,2)  EN(I),(RN(J,I),Js1,9) 

2  FORMAT (6X,F 9. 3,9 (3X,E9.3 ) ) 

3  CONTINUE 

IF  ENERGY  FUG  SET,  RETURN  TO  ENERGY  PART. 

IF(NTRY.EQ. 1 )  GO  TO  100 
5  CONTINUE 

IZsZ+,5 


FIND  THE  NEAREST  TABUUTED  ION. 

KsKK(IZ) 

COMPUTE  THE  APPROXIMATE  INDEX  IN  THE  TABLE. 
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IJ  =66 

IF(E.GT. 8000, )  GO  TO  10 
IJs(AL0G (E )+4.605 ) *4. 85567-1. 0 
IJaMAXOdJ,  D 
IJ  =MIN0 (I J,66  ) 

FIND  THE  TABULATED  ENERGY  JUST  ABOVE  E. 

10  CONTINUE 

DO  11  I=IJ,66 
IF(EN(I).LT.E)  GO  TO  11 
J=I 

GO  TO  13 

11  CONTINUE 

EXTRAPOLATE  THE  RANGE  ACCORDING  TO  THE  PLATEAU 
VALUE  OF  THE  STOPPING  POWER. 

RNK=RN (K ,  66 )+  (E-EN  (66  ) )  /DEDX(K  ) 

GO  TO  14 

COMPUTE  THE  INTERPOLATION  FACTOR. 

13  CONTINUE 

F  =  (E-EN(J-1  ))/(EN(J  )-EN(J-1  )) 

FsAMINKF,  1.0) 

INTERPOLATE  RANGE  IN  THE  TABLE. 

RNKsF* (R N  (K ,  J  )-RN (K ,  J-1  ))+RN(K,  J-1  ) 

INTERPOLATE  TO  THE  ION  IN  QUESTION. 

14  CONTINUE 

RAL=(A/( Z*Z) )* (ZN (K )*ZN(K )/AN(K ) )*RNK 
RETURN 

ENTRY  EAL(Z,A,R) 

THIS  ENTRY  RETURNS  THE  ENERGY  (IN  MEV/AMU)  OF  AN  ION  HAVING 
A  RESIDUAL  RANGE  R,  IN  G/CM««SQ. 

Z  IS  THE  ION  ATOMIC  NUMBER. 

A  IS  THE  ION  ATOMIC  MASS. 

SET  ENTRY  FLAG  FOR  ENERGY. 

IZ=Z+.  5 
NTRY=1 

IF  THIS  IS  THE  FIRST  CALL,  READ  IN  THE  TABLES. 

IF(IST.EQ.O)  GO  TO  1 
C 
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C  FIND  THE  NEAREST  TABULATED  ION. 

C 

100  CONTINUE 
KsKK(IZ) 

COMPUTE  THE  ION  INTERPOLATION  FACTOR. 
110  CONTINUE 

FK  =  (A/(Z*Z))*(ZN(K)*ZN(K)/AN  (K)> 

FIND  THE  TABULATED  RANGE  JUST  ABOVE  R. 

DO  120  1=1,66 
RT=FK«RN(K,I) 

IFCRT.LT. R)  GO  TO  120 
J=I 

GO  TO  130 
120  CONTINUE 

EXTRAPOLATE  THE  ENERGY. 

RAT=DEDX(K )/FK 

EAL=EN  (66)+RAT*(R-FK#RN  (K,66  ) ) 

RETURN 

INTERPOLATE  THE  ENERGY. 

130  CONTINUE 

F  =  (R-FK*RN  (K,  J-1  )  )/(RT-FK*RN(K,  J-1  )) 
EAL  =F*  ( EN  ( J  )  -EN  ( J  -1  ) )  * EN  ( J  -1  ) 

RETURN 

END 
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FUNCTION  CUT  (Z,E,Y,M) 

THIS  ROUTINE  OBTAINS  DIFFERENTIAL  PARTICLE  FLUXES  FROM  CRF, 
APPLIES  THE  GEOMAGNETIC  CUTOFF  TRANSMISSION  FUNCTION 
IN  THE  GTRANS.DAT  FILE,  AND  RETURNS 

THE  RESULTING  FLUX,  MODULATED  TO  THE  ORBIT -AVERAGE  CUTOFF. 

THIS  ROUTINE  ALSO  ADDS  SINGLY-IONIZED  ANOMALOUS  COMPONENT  NUCLEI 
TO  THE  COSMIC  RAY  SPECTRA  (WHEN  M=4). 

Z  =  ION  ATOMIC  NUMBER. 

E  =  ION  ENERGY  IN  MEV/AMU. 

Y  =  YEAR  (1975.144  =  SOLAR  MIN. j  1980.598  *  SOLAR  MAX.). 

M  s  INTERPLANETARY  MEDIUM  WEATHER  INDEX. 

COMMON /MASS/ A( 92) 

DIMENSION  R (200) ,T (200  ) 

DATA  ITS/-1 / 

IF(ITS.EQ.  1  )  GO  TO  6 

OPEN (UNIT=1 5 , READONLY, SHARED, STATUS  =  'OLD ' ,FILE  =  'GTRANS. DAT ' ) 

ITS  si 

DO  5  I  si, 200 

READ (1 5, 1  )  R  (I)  ,T  (I ) 

FORMAT (5X,F6. 3, 5X,F8. 6) 

CONTINUE 
SHADOW  =T  (200) 

IZ  =Z+.  5 
CUTsO.O 
ESsE 


DON'T  CALL  CRF  WITH  ENERGIES  BELOW  10  MEV/U. 

IF(E.LT. 10.)  ES=10. 

MK=M 

MQ=M 


Ms4  MEANS  A  SINGLY-IONIZED  ANOMALOUS  COMPONENT  +  GALACTIC  COSMIC 
RAYS. 

GET  THE  GALACTIC  COSMIC  RAYS  HERE,  THE  ANOMALOUS  COMPONENT  WILL 
BE  ADDED  BELOW. 

IF(M.EQ.  4)  MQsl 

GET  THE  DIFFERENTIAL  FLUX. 

FsCRF(Z,ES,Y,MQ) 

COMPUTE  THE  MAGNETIC  RIGIDITY. 

Ps(A(IZ)  /Z)  *(ES*ES+1 862.  324»ES)««  .5/1 000 

LOOK  UP  THE  TABULATED  MAGNETIC  RIGIDITY  JUST  ABOVE  P. 


4  DO  2  1=2,200 

IFCP.GT. R(I  ))GOTO  2 
ISAV  =1 
GOTO  3 

2  CONTINUE 
TRsSHADOW 
GO  TO  10 

C 

C  INTERPOLATE  THE  TRANSMISSION  FACTOR  (AVERAGED  FOR  THE  ORBIT). 

C 

3  TR  =  (T(ISAV)-T  (ISAV-1  ))» (P-R  (ISAV-1  )) /(R  (ISAV)-R  (ISAV-1  )) 

1  +T(  ISAV-1) 

C 

C  APPLY  THE  TRANSMISSION  FACTOR  TO  MODULATE  THE  FLUX. 

C 

10  CUT=CUT+F*TR 

C 

C  ADD  IN  THE  DIRECT  CONTRIBUTION  FROM  THE  TRAPPED  PROTONS. 

C 

IF(IZ.EQ.I)  CUT=CUT+PROTON  (ES ) 

C 

C  UNLESS  A  SINGLY-IONIZED  ANOMALOUS  COMPONENT  MUST  BE  ADDED,  RETURN. 

C 

IFCMK.NE.4)  GO  TO  12 
C 

C  ADD  SINGLY-IONIZED  HELIUM. 

C 

IFCIZ.NE.2)  GO  TO  40 
F=(1.54E+4)«ES»*(-2) 

IF(ES. LT. 195. )  F=. 4 
P=4.»(ES«ES+1 862. 324  «ES)*«. 5/1000 
MK=1 
GO  TO  4 
C 

C  ADD  SINGLY-IONIZED  CARBON. 

C 

40  IF(IZ.NE.6)  GO  TO  45 

IF(ES. LT. 10.)  GO  TO  43 
F  =0. 27*ES**(-2 ) 

GO  TO  44 

43  F  =  (4 .  OE-3  )*EXP  (-  (ALOG  (ES  )-1 . 79  )**2  /. 7 ) 

IF(F. LT. 0. )  F=0 

44  P=12!«(ES«ES+1862.324*ES)»».5/1000. 

MK  =  1 

GO  TO  4 
C 

C  ADD  SINGLY-IONIZED  NITROGEN. 

C 

45  IFQZ.NE.7)  GO  TO  50 
IF (ES. LT.20. )  GO  TO  41 
F=.773*ES»*(-2) 

GO  TO  42 


F=(1.54E-2)*EXP(-(ALOG(ES)-1.79)**2/.7) 
IF ( F.  LT.  0. )  FsO. 

P=14.»(ES*ES+1862«324  *ES)*». 5/1000. 

MK=1 
GO  TO  4 

ADD  SINGLY-IONIZED  OXYGEN. 

IF(IZ.NE.8)  GO  TO  60 
IF(ES.LT.30.)  GO  TO  51 
F=1 . 32*ES**(-2. ) 

GO  TO  52 

F  =  (6.  OE-2  )*EXP  (- (ALOG(ES  )-1 . 79  )**2/.7) 
IFCF.LT.O.)  F =0, 

Psl  6*<ES«ES+1  862.  324 «ES  )*»  .5/1 000. 

MK  =  1 
GO  TO  4 

ADD  SINGLY-IONIZED  NEON. 

IFUZ.NE.10.)  GO  TO  70 
IF(ES.LT.20. )  GO  TO  61 
F=0.4«ES«»(-2) 

GO  TO  62 

F  =  (8 .  OE-3  )*EXP  (-  ( ALOG  (ES  )-1 . 79  )*»2  / . 7 ) 
IF(F.LT. 0.)  F  =0 . 

P=20*(ES*ES+1  862.  324*ES)**.5/1 000. 

MK  =  1 
GO  TO  4 

ADD  SINGLY-IONIZED  MAGNESIUM. 

IFCIZ.NE.12)  GO  TO  80 
IF(ES. LT.20. )  GO  TO  71 
F  =0. 16#ES*#(-2 ) 

GO  TO  72 

F  =  (8. OE-4  )*EXP  (-(ALOGCES  )-2.  3)**2/.7 ) 
IFCF.LT.O.)  FsO. 

Ps24. 3*(ES«ES+1 862. 32^*ES)*».5/1000. 

MK  =  1 
GO  TO  4 


ADD  SINGLY-IONIZED  SILICON. 

IF(IZ.NE.14)  GO  TO  85 
IF(ES. LT. 10. )  GO  TO  81 
FsO,  1  *ES**(-2 ) 

GO  TO  82 

F  =  (1 .  OE-3  )*EXP  C-(AL0G(ES)-2.2)«*2/0.4) 
IF(F. LT. 0.)  FsO. 

Ps28. * (ES*ES+1 862.  324  «ES  )**  .5/1  000. 

MK  si 
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GO  TO  4 


ADD  SINGLY-IONIZED  ARGON. 

85  IFCIZ.NE. 18)  GO  TO  90 
IF(ES, LT.20. )  GO  TO  86 
F=0. 028*ES**C-2 ) 

GO  TO  87 

86  F  =  (5.4E-4)*EXP(-(AL0G(ES)-1 .79  )**2/.7) 

IF( F. LT  0  )  F  =0 

87  P=40*.«(ES«ES+1  862.  324«ES)*». 5/1000. 

MK  =  1 

GO  TO  4 

ADD  SINGLY-IONIZED  IRON. 

90  IF(IZ.NE.26)  GO  TO  12 

IFCES.LT. 30. )  GO  TO  91 
F=0. 35*ES**(-2 ) 

GO  TO  92 

F  =  (6 . OE-4 )*EXP  C-CALOGCES )-2.  48  )**2/2. ) 

IFCF.LT.O.)  F =0. 

P=56. #(ES*ES+1 862. 324 *ES )*• .5/1 000. 

MK  =  1 
GO  TO  4 

THERE  PROBABLY  ARE  ANOMALOUS  COMPONENTS  IN  THE  SPECTRA  OF 
THE  NUCLEI  HEAVIER  THAN  IRON,  BUT  THERE  ARE  NO  DATA  ON  THEM 
AT  THIS  TIME,  SO  THEY  CANNOT  BE  INCLUDED  IN  THE  MODEL. 

12  RETURN 
END 
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PROGRAM  BENDEL 


THIS  PROGRAM  USES  BENDEL 'S  FORMULA  ("Proton  Upsets  in  Orbit",  by 
W.  L.  Bendel  and  E.  L.  Petersen,  IEEE  Trans  in  Nucl.  Sci.,  Vol.  NS-30, 
p.  4481-5,  1983)  TO  CALCULATE  THE  UPSET  RATE  CAUSED  BY  PROTONS.  THE 
PARAMETER  A  MUST  BE  INPUT.  THIS  PROGRAM  WAS  CREATED  BY  MODIFYING  THE 
SPEC  PROGRAM.  SUBROUTINESi  INSIDE,  CUT,  CRF,  DEDXAL,  RAL. 

REAL  INSIDE 
Z  =  1. 

WRITE  (6,103) 

103  FORMAT ('  ENTER  THE  NUMBER  OF  ENERGY  STEPS  USED  IN  THE'/ 

1'  CALCULATION  (MAXIMUM:  1000  ):  ',$) 

ACCEPT  *,NDIF 
WRITE  (6,110) 

110  FORMAT ( '  ENTER  THE  YEAR  FOR  WHICH  YOU  WANT  THE  PROTON -INDUCED  »/ 

1'  UPSET  RATE  (1975.144  FOR  SOLAR  MIN.;  1980.598  FOR  SOLAR  MAX.)', 

2':  ',$) 

ACCEPT  «,Y 
WRITE  (6,111) 

111  FORMAT ('  ENTER  THE  ORBITAL  INDEX:  0  FOR  NO  CUTOFF  AND  NO  ’, 

1  'TRAPPED  PROTONS;'/'  1  FOR  CUTOFF  AND  TRAPPED  PROTONS', 

2'  (NOTE:  STASS.DAT  AND  GTRANS.DAT'/ 

3'  MUST  BE  SUPPLIED  WHEN  1  IS  ENTERED  HERE):  ',$) 

ACCEPT  *,IT 
WRITE  (6,112) 

112  FORMAT (  '  ENTER  THE  INTERPLANETARY  WEATHER  INDEX:  ',$) 

ACCEPT  * ,M 

WRITE  (6,113) 

113  FORMAT ('  ENTER  THE  SHIELDING  THICKNESS  IN  INCHES  OF  »/ 

1'  ALUMINUM  (OR  EQUIVALENT):  ',$> 

ACCEPT  *,THK 
WRITE  (6,114) 

114  FORMAT ( '  ENTER  THE  PARAMETER  A  (ME7) .  (ENTER  0'/ 

1'  IF  YOU  DON '  'T  KNOW  IT):  ',$) 

ACCEPT  «,A 

IF  (A.GT.O.)  GOTO  150 

IF  NOTHING  WAS  ENTERED,  FIGURE  OUT  THE  VALUE  OF  ''A''  FROM 
EXPERIMENTAL  DATA. 

WRITE  (6,115) 

115  FORMAT ('  ENTER  AN  EXPERIMENTALLY  MEASURED  UPSET  CROSS-SECTION'/ 

1'  IN  UNITS  OF  (UPSETS/BIT )/(PR0T0N/CM**2 ):  ',$) 

ACCEPT  *,XSECT 
WRITE  (6,116) 

116  FORMAT ('  ENTER  THE  ENERGY  (MEV)  AT  WHICH  THIS  CROSS-SECTION'/ 

1'  MEASUREMENT  WAS  MADE:  ',$) 

ACCEPT  *,  ENERGY 
ALOW  =  1. 

A HIGH  =  200. 
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XSECT  *  1,E12*XSECT 
DO  130  J=1 , 16 

A  =  .5*(AHIGH+AL0W) 

X  =  SIGMA( A, ENERGY) 

SIGMA  IS  A  DECREASING  FUNCTION  OF  A. 

IF  (X.GT. XSECT)  ALOW=A 
IF  (X.LE. XSECT)  AHIGH=A 
130  CONTINUE 

150  CONTINUE 

CONVERT  THE  THICKNESS  TO  G/CM«*2  (ALUMINUM  EQUIVALENT). 

TH=THK*6. 8528 

TYPE  OUT  THE  INPUTTED  DATA. 

WRITE  (6,117)  Y 

117  F0RMAT(//'  THE  UPSET  RATE  WILL  BE  FOR  THE  YEAR  ,,F8.3,,.f) 

IF(IT.EQ.O)  WRITE  (6,118) 

118  FORMAT ('  THE  COSMIC  RAYS  WILL  BE  UNSHIELDED  BY  THE',/, 

1  1  EARTH ' 'S  MAGNETIC  FIELD,  AND  THE  CONTRIBUTION  FROM  TRAPPED',/, 

2  '  PROTONS  WILL  NOT  BE  INCLUDED. ' ) 

IF(IT.EQ. 1 )  WRITE  (6,119) 

119  FORMAT ('  THE  UPSET  RATE  WILL  INCLUDE  THE  EFFECTS  OF',/, 

1  *  GEOMAGNETIC  SHIELDING  ACCORDING  TO  THE  GEOMAGNETIC  CUTOFF',/, 

2  '  TRANSMITTANCE  FUNCTION  TABULATED  IN  GTRANS.DAT.  THE',/, 

3  '  CONTRIBUTION  FROM  TRAPPED  PROTONS  WILL  BE  INCLUDED',/, 

H  '  ACCORDING  TO  THE  ORBIT-AVERAGED  TRAPPED  PROTON  SPECTRA',/, 

5  '  TABULATED  IN  STASS.DAT.') 

WRITE  (6,120)  M 

120  FORMAT ('  THE  INTERPLANETARY  WEATHER  INDEX  IS  ',12,'.') 

WRITE  (6,121 )  TH 

121  FORMAT ('  THE  SHIELDING  THICKNESS  IS  ',F7.3,'  GM/CM«»2 
1  ALUMINUM  (OR  EQUIVALENT).') 

WRITE  (6,122)  A 

122  FORMAT ( '  THE  UPSET  THRESHOLD  PARAMETER  IS  ',F6.2,'  MEV. ' ) 

SET  UP  THE  ENERGY  RANGE  AND  DIVIDE  IT  INTO  EQUAL  LOGARITHMIC  STEPS. 

EMAX=32000. 

EMINslO. 

ELNX=ALOG (EMIN) 

DELN = ( ALOG  ( EMAX )  -ELNX )  /ND I F 
UPSETS =0. 

COMPUTE  THE  DIFFERENTIAL  FLUX  AT  EACH  ENERGY.  THE  SUBPROGRAMS  RETURN 
FLUX  IN  PARTICLES/((M*»2)«STER»SEC)/(MEV/(G/CM«»2)). 

DO  200  Isl.NDIF 
EN sEXP (ELNX) 
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«.*'.  kj *^.  ™*r*  «  v  .-w  *■_■»■■  -*  p.«v  w™  -»<»-■»»  .-u.:  -,Jl™.nUiJ-w 


ELNX=ELNX+DELN 
F  =INSIDE (Z,EN, Y,M,IT,TH ) 

COMPUTE  BENDEL*S  FLUX.  E  AND  A  ARE  IN  MEV,  Y  IS  DIMENSIONLESS, 
SIGMA  IS  IN  IE-12  (UPSETS/BIT )/(PROTON/CM**2) ,  UPSETS  IS  IN 
UPSETS/ BIT»DAY. 

THE  FACTORS  IN  THE  CONVERSION  CONSTANT  ARE: 

IE-12  FROM  THE  DEFINITION  OF  SIGMA; 

86400  TO  CONVERT  FROM  SECONDS  TO  DAYS; 

.01 **2  TO  CONVERT  FROM  M*»2  TO  CM*«2; 

4*PI  TO  GET  RID  OF  THE  STE RADIANS. 

CONST =1 . E— 1 2  * ( .01 **2. ) *(4 . *3. 1416) 

UPSETS  =  UPSETS  +  CONST »SIGMA(.A,EN)*F«  (EXP  (ELNX)-EN  ) 

200  CONTINUE 

WRITE (6,14  )UPSETS 

14  FORMAT ('  ERROR  RATE  =  ',1PE12.5,'  IN  UPSETS/EIT*SEC. ‘ ) 

UPSETS  =UPSETS  *86400 . 

WRITE (6, 15 )UPSETS 

15  FORMAT (’  ERROR  RATE  =  ,,1PE12.5,*  IN  UPSETS/ BIT«DA Y. * ) 

END 
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FUNCTION  SIGMA  ( A,E ) 

THIS  IS  THE  SIGMA  FROM  PAGE  4484  OF  THE  BENDEL  PAPER. 

Y  =  (SORT  (18./A))*(E-A) 

IF  (Y.LT.O.)  Y=0. 

SIGMA  =  ((24./A)**l4.)»((1.-EXP(-.l8*SQRT(Y))}»»4.) 

RETURN 

END 
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PROGRAM  SFEC 

THIS  PROGRAM  USES  THE  SUBPROGRAMS  INSIDE,  CUT,  CRF,  DEDXAL, 

AND  RAL  TO  TABULATE  THE  DIFFERENTIAL  AND  INTEGRAL  ENERGY  SPECTRA  OF 
COSMIC  RAYS  OF  A  SPECIFIED  CHARGE. 

REAL  INSIDE 

COMMON  /MASS/AMSTAE(92) 

CHARACTER* 12  ALPHA, BETA, GAMMA 
DIMENSION  F< 1000) ,E (1001 ) 

DATA  ND IF/0 /.NINT/O/, GAMMA/*  */ 

WRITE  (6,101) 

101  FORMAT ( *  ENTER  THE  NAME  OF  THE  DIFFERENTIAL  ENERGY  SPECTRUM',/, 

1*  (ENTER  NOTHING  IF  YOU  DO  NOT  WANT  ONE):  »,$) 

ACCEPT  102, ALPHA 

102  F0RMAT(A12 ) 

IF ( ALPHA. EQ. GAMMA)  GO  TO  104 
OPEN (UNIT=40,FILE=ALPHA ) 

WRITE  (6,103) 

103  FORMAT ('  ENTER  THE  NUMBER  OF  POINTS  IN  THE  TABULATION  OF',/, 

1'  THE  DIFFERENTIAL  ENERGY  SPECTRUM  (MAXIMUM:  1000):  ',$) 

ACCEPT  *,NDIF 

104  WRITE  (6,105) 

105  FORMAT ( '  ENTER  THE  NAME  OF  THE  INTEGRAL  ENERGY  SPECTRUM',/, 

1'  (ENTER  NOTHING  IF  YOU  DO  NOT  WANT  ONE):  ',$) 

ACCEPT  102, BETA 

IF( BETA. EQ. GAMMA)  GO  TO  107 

OPEN  (UNITs45 ,  FILE=BETA  ) 

WRITE  (6,106) 

106  FORMAT ('  ENTER  THE  NUMBER  OF  POINTS  IN  THE  TABULATION  OF',/, 

1  '  THE  INTEGRAL  ENERGY  SPECTRUM  (MAXIMUM:  1000):  ',$) 

ACCEPT  *,NINT 

107  NPT  =MAX0 (NDIF, NINT ) 

WRITE  (6,108) 

108  FORMAT ('  ENTER  THE  ATOMIC  NUMBER  OF  THE  ELEMENT  FOR  WHICH’,/, 

1*  YOU  WANT  THE  ENERGY  SPECTRUM:  ',$) 

ACCEPT  *,IELM 
Z=IELM 

WRITE  (6,110) 

110  FORMAT ( '  ENTER  THE  YEAR  FOR  WHICH  YOU  WANT  THE  SPECTRUM',/, 

1'  (1975.144  FOR  SOLAR  MIN.;  1980.598  FOR  SOLAR  MAX..):  ',$) 

ACCEPT  *,Y 
WRITE  (6,111) 

111  FORMAT ('  ENTER  THE  ORBITAL  INDEX:  0  FOR  NO  CUTOFF  AND  NO  ', 

1  'TRAPPED  PROTONS?'/'  1  FOR  CUTOFF  AND  TRAPPED  PROTONS', 

2'  (NOTE:  STASS.DAT  AND  GTRANS. DAT ' / 

3'  MUST  BE  SUPPLIED  WHEN  1  IS  ENTERED  HERE):  ',$) 

ACCEPT  *,IT 
WRITE  (6,112) 

112  FORMATC  ENTER  THE  INTERPLANETARY  WEATHER  INDEX:  ',$) 

ACCEPT  * ,M 
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WRITE  (6,113) 

113  FORMAT ('  ENTER  THE  SHIELDING  THICKNESS  IN  INCHES  OF  ’,/, 

1*  ALUMINUM  (OR  EQUIVALENT):  ',$) 

ACCEPT  *,THK 

CONVERT  THE  THICKNESS  TO  G/CM**2  (ALUMINUM  EQUIVALENT). 

TH=THK#6 , 8528 
WRITE  (6,116)  IELM 

116  FORMAT (//'  THE  ENERGY  SPECTRA  WILL  INCLUDE  THE  CONTRIBUTION  1 / 

1  '  OF  ELEMENT  ',12,'  .') 

IF(NDIF.GT.O)  WRITE  (6,114)  ALPHA, NPT 
IF(NINT.GT.O)  WRITE  (6,115)  BETA, NPT 

114  FORMAT ('  THE  FILE  ',A12,'  WILL  CONTAIN  THE  DIFFERENTIAL',/ 

1  ,'  ENERGY  SPECTRUM,  TABULATED  IN  i,l4,»  DATA  POINTS.*) 

115  FORMAT ( '  THE  FILE  ',A12,'  WILL  CONTAIN  THE  INTEGRAL',/ 

1  ,'  ENERGY  SPECTRUM,  TABULATED  IN  ',14,'  DATA  POINTS.’) 

WRITE  (6,117)  Y 

117  FORMAT ('  THESE  ENERGY  SPECTRA  WILL  BE  FOR  THE  YEAR  ’,F8.3,'.') 
IF(IT.EQ.O)  WRITE  (6,118) 

118  FORMAT ('  THESE  ENERGY  SPECTRA  WILL  BE  UNSHIELDED  BY  THE',/, 

1  '  EARTH ' 'S  MAGNETIC  FIELD,  AND  THE  CONTRIBUTION  FROM  TRAPPED',/, 

2  '  PROTONS  WILL  NOT  BE  INCLUDED. ' ) 

IF(IT.EQ.I)  WRITE  (6,119) 

119  FORMAT ('  THESE  ENERGY  SPECTRA  WILL  INCLUDE  THE  EFFECTS  OF’,/, 

1  '  GEOMAGNETIC  SHIELDING  ACCORDING  TO  THE  GEOMAGNETIC  CUTOFF',/, 

2  '  TRANSMITTANCE  FUNCTION  TABULATED  IN  GTRANS.DAT.  THE’,/, 

3  '  CONTRIBUTION  FROM  TRAPPED  PROTONS  WILL  BE  INCLUDED',/, 

4  '  ACCORDING  TO  THE  ORBIT-AVERAGED  TRAPPED  PROTON  SPECTRA’,/, 

5  '  TABULATED  IN  STASS.DAT. ') 

WRITE  (6,120)  M 

120  FORMAT (*  THE  INTERPLANETARY  WEATHER  INDEX  IS  ',12,'.') 

WRITE  (6,121 )  TH 

FORMAT ('  THE  SHIELDING  THICKNESS  IS  »,F7.3, '  GM/CM«»2 
1  ALUMINUM  (OR  EQUIVALENT).',//) 

SET  UP  THE  ENERGY  RANGE  AND  DIVIDE  IT  INTO  EQUAL  LOGARITHMIC  STEPS. 

EMAX=32000. 

EMIN  =  10. 

ELNX=ALOG (EM IN ) 

DELN  s  ( ALOG  ( EMAX )  -ELNX )  /N  PT 

COMPUTE  THE  DIFFERENTIAL  FLUX  AT  EACH  ENERGY. 

DO  200  I  si ,  NPT 
EN  =EXP(ELNX) 

ELNX=ELNX+DELN 
E (I )sEN 

F(I)=INSIDE(Z,EN,Y,M,IT,TH) 

IF(NDIF.GT.O)  WRITE  (40, 201 )  F(I),E(I) 

FORMAT (2X, E 12. 5, 2X, E12,  5 ) 


201 


CONTINUE 


SUM  THE  DIFFERENTIAL  ENERGY  SPECTRUM  TO  GET  THE  INTEGRAL  ENERGY 
SPECTRUM  (BOX  INTEGRATION) , 

SUM=. 6666*1 NSIDE (Z,EN, Y,M, IT, TH  )*EN 
DO  300  I=NPT,2,-1 

SUMsSUM+0. 5*(F(I  )+F(I-1  ) )*(E  (I  )-E  (1-1  )) 

IF(NINT.GT.O)  WRITE  (45, 201  )  SUM.ElI-l) 

CONTINUE 
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PROGRAM  UPSET 

C  THIS  PROGRAM  COMPUTES  THE  UPSET  RATE  DUE  TO  THE 

C  DIRECT  IONIZATION  OF  INDIVIDUAL  PARTICLES.  IT  ASSUMES 

C  THAT  FOR  EACH  BIT  THERE  EXISTS  A  SENSITIVE  VOLUME. 

C  IF  AN  AMOUNT  OF  ELECTRICAL  CHARGE  OQCRIT)  IS  CREATED 

C  WITHIN  THIS  VOLUME  BY  THE  IONIZING  PARTICLE,  THEN  AN 

C  UPSET  WILL  OCCUR.  THE  SENSITIVE  VOLUME  IS  IDEALIZED  AS 

C  A  PARALLELEPIPED  WITH  DIMENSIONS  X,  Y,  AND  Z. 

C  THE  INPUTS  ARE: 

C  THE  CRITICAL  CHARGE  QCRIT  IN  PICOCOULOMBS, 

C  THE  DIMENSIONS  OF  THE  SENSITIVE  VOLUME  IN  MICROMETERS, 

C  THE  NAME  OF  THE  FILE  CONTAINING  THE  INTEGRAL  LET 

C  SPECTRUM  (SPEC. INT  IS  ASSUMED  IF  NO  NAME  IS  ENTERED). 

C  THE  INTEGRAL  LET  SPECTRUM  IS  ASSUMED  TO  BE  IN  UNITS  OF 

C  PARTICLES/ (SQUARE  METER *STERAD IAN *SEC  )  VERSUS 

C  MEV* (SQUARE  CENTIMETER) /GRAM. 

C  THE  OUTPUT  IS  GIVEN  IN  UPSETS/( BIT "SECOND )  AND  UPSETS/(BIT*DAY) . 

THIS  CALCULATION  USES  THE  METHOD  DESCRIBED  IN  "THE  VARIABILITY 
OF  SINGLE  EVENT  UPSET  RATES  IN  THE  NATURAL  ENVIRONMENT," 

JAMES  H.  ADAMS,  JR.,  IEEE  TRANS.  ON  NUCL.  SCI.,  NS-30,  4475- 
4480,  DEC.,  1983. 

IT  IS  OPERATIONALLY  EQUIVALENT  TO  ROCKWELL’S  CRIER  PROGRAM 
(PICKEL  AND  BLANDFORD,  IEEE  TRANS.  ON  NUCL.  SCI.  NS-26,  DEC. 
1979,  PP. 4735-4739)  WHEN  USED  WITH  HEINRICH’S  LET  SPECTRUM 
(W.  HEINRICH,  RADIATION  EFFECTS,  VOL.  34,  PP.  143-8,  1977). 

THIS  PROGRAM  CALLS  A  SUBPROGRAM,  DIFPLD,  THAT  RETURNS 
THE  DIFFERENTAL  PATHLENGTH  DISTRIBUTION  IN  THE  SENSTIVE 
VOLUME. 

CHARACTERS  ALPHA, DEFAULT, BLANK 
REAL  L  (1000)  ,LMIN 
DIMENSION  FLUX (1000) 

DATA  DSI/2. 321 /.BLANK/’  ' /.DEFAULT/ ’SPEC. INT'/ 

GET  THE  INPUT  DATA. 

WRITE (6, 1) 

FORMAT ('  ENTER  THE  CRITICAL  CHARGE  IN  PICOCOULOMBS:  ’,$) 

ACCEPT  *,  QCRIT 
WRITE (6, 2) 

FORMAT ('  ENTER  THE  DIMENSIONS  OF  THE  SENSITIVE  VOLUME, 

1  X,  Y,  Z,  IN  MICROMETERS:  ',$) 

ACCEPT  »,  X.Y.Z 

CONVERT  FROM  MICROMETERS  TO  CENTIMETERS. 

X=X».0001 
YsY«.0001 
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Z=Z*.0001 
WRITE(6, 3) 

3  FORMAT (•  ENTER  THE  NAME  OF  THE  INTEGRAL  LET  FILE:  ',$) 
ACCEPT  4,  ALPHA 

4  FORMAT (A  12) 

IF  NO  FILE  NAME  WAS  ENTERED,  ASSUME  'SPEC.INT*. 

IF ( ALPHA. EQ. BLANK)  ALPHA=DEFAULT 

OPEN (UNIT  =  1 , READONLY, SHARED, STATUS s 'OLD ' .FILE  sALPHA ) 

READ  THE  INTEGRAL  LET  SPECTRUM. 

FLUX  IS  IN  PARTICLES/ (SQUARE  METER »STERAD IAN *SEC ) . 

L  IS  IN  MEV* (SQUARE  CENTIMETER )/GRAM. 

N=1 

5  CONTINUE 

READ  (1 , 1 1,  END«6  )  FLUX(N),L(N) 

11  FORMAT  (2X,E12.5,2X,E12.5 ) 

N=N+1 

DON'T  OVERFILL  THE  ARRAYS. 

IF(N.LE. 1000)  GOTO  5 

COMPUTE  THE  NUMBER  OF  POINTS  IN  THE  LET  SPECTRUM. 

6  CONTINUE 
NsN-1 

COMPUTE  THE  SURFACE  AREA  OF  THE  SENSITIVE  VOLUME. 

Ss  (2 .  *X  *Y+2 .  *X  *Z+2 .  *Y  *Z) 

CONVERT  FROM  SQUARE  CENTIMETERS  TO  SQUARE  METERS. 

S=S».0001 

CONVERT  THE  DIMENSIONS  OF  THE  SENSITIVE  VOLUME  TO  G/CM**2. 

XsX«DSI 

YsY»DSI 

Z=Z*DSI 

COMPUTE  THE  MAJOR  DIAMETER  OF  THE  SENSITIVE  VOLUME. 

PMAXsSQRT  (X  *X+Y  »Y+Z»  Z) 

COMPUTE  THE  ENERGY  (IN  MEV)  REQUIRED  TO  PRODUCE  QCRIT(IN  PC) 
HOLE-ELECTRON  PAIRS  IN  SILICON. 


ENERGY s22. 5 *QCRIT 
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L\P 

8 


u 


c 

C  COMPUTE  THE  MINIMUM  LET  THAT  CAN  PRODUCE  AN  UPSET. 

C 

LMINsENERGY/PMAX 
SUM =0.0 
QsLMIN 
C 

C  INTEGRATE  FROM  LM IN  TO  THE  LARGEST  LET  IN  THE  SPECTRUM. 

C 

DO  10  I  *N ,  1 ,  -1 
IF(L (I).LT.LMIN )  GO  TO  10 
C 

C  COMPUTE  THE  PATHLENGTH  CORRESPONDING  TO  L(I). 

C 

D=ENERGY/L  (I ) 

C 

C  CARRY  OUT  THE  INTEGRAL. 

C 

SUM=SUM+(L  (I  )-Q  )«DIFPLD  (D ,  X,  Y, Z)  »FLUX(I  )/(L  (I  )«*2 ) 
Q=L(I) 

10  CONTINUE 

C 

C  COMPUTE  THE  ERROR  RATE. 

C 

12  CONTINUE 

ERR=ENERGY«S*3.1416»SUM 
WRITE (6, 14 )ERR 

14  FORMAT ('  ERROR  RATE  =  ’.IPEte.S,’  IN  UPSETS/ BIT «SEC. ' ) 

C 

C  CONVERT  TO  UPSET S/ BIT »DAY. 

C 

ERR=ERR«86400. 

WRITE (6, 15 )ERR 

15  FORMAT ( '  ERROR  RATE  =  ’.IPEte.S,’  IN  UPSETS/ BIT  «DA  Y. 1 ) 
END 
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FUNCTION  DIFPLD(S,L,W,H) 

THIS  FUNCTION  RETURNS  THE  PROBABILITY  DENSITY  FOR  PATHS 
OF  LENGTH  S  THROUGH  A  PARALLELEPIPED  OF  DIMENSIONS 
L,  W,  AND  H.  S,  L,  W,  AND  H  MUST  BE  IN  THE  SAME  UNITS. 

THIS  IS  AN  EXACT  SOLUTION,  DUE  TO  M.  D.  PETROFF  OF 
ROCKWELL  INTERNATIONAL  (SEE  J.  C.  PICKEL  AND  J.  T.  BLANDFORD, 
IEEE  TRANS.  ON  NUCL.  SCI.  NS-27,  1006(1980))  WITH 
SIMPLIFICATIONS  DUE  TO  WARREN  BENDEL  OF  NRL  (PRIVATE 
COMMUNICATION).  THE  EQUATION  NUMBERS  REFER  TO  THE  APPENDIX 
OF  PICKEL  AND  BLANDFORD'S  PAPER. 

REAL  L 

EQUATION  (A-7) 

•» 

AP*3.#(H«W4H«L+L»W) 

EQUATION  (A-8) 

DIFPLDs (G (S , L,W , H )+G  (S,W,L,H)+G  (S,L,H,W)+G(S,W,H,L)+ 

1  G(S,H,W,L)4G(S,H,L,W))/(3.1J»16»AP) 

RETURN 

END 
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FUNCTION  G(S,Xt Y-Z) 
REAL  KSQ 


PRELIMINARY  DEFINITIONS 

KSQ=X*X+Y*Y 
TSQsX*X+Z*Z 
TsSQRT  (TSQ) 

RSQ*KSQ+Z»Z 

RsSQRT(RSQ) 

Vs12.*X*Y*Z*Z 
PSQ=S«S-Z«Z 
QSQ  sS  *S-X  *X  -Z#  Z 

IF((S.GE.0.0).AND. (S.LT.Z))  GO  TO  10 
IF(  (S.GE.Z)  .AND.  (S.  LT»T )  )GO  TO.  20 
IF ( ( S . GE . T )  .AND. (S. LE. R ) )G0  TO  30 
G=0.0  ... 

RETURN 

EQUATION  (A-9) 

10  G=8.*Y*Y*Z/KSQ-S*(3.*X*Y/(R*T  ))**2 

RETURN 


EQUATION  (A-10) 

G=S*(3.*Y/SQRT  (KSQ)  )**2-S*(3.*X*Y/(T*R  )  )**2 

1  -X*  (SQRT  (PSQ)/S )*(8.  +4.#Z*Z/  (S*S ) ) 

2  + (V  « ATA  N  (Y  /X  )-  (Y  *Z«  Z/ SQRT  (KSQ  ) ) »«2 )  /  (S  «S  «S  ) 
RETURN 

EQUATION  (A-11) 

30  Gs~S*(3.*X*Z/(R*SQRT  (KSQ))  )**2 

1  +  (X  «X*  1*  l*  ( Z«  Z/KSQ-3 .  )+V*ATAN  (Y  /X ) )  /  (S  «S  «S  ) 

2  +Y*Z*  Z*(SQRT  (QSQ)/S)*(8./TSQ+4,/(S*S ) ) 

3  -(V/(S*S*S)  )*ACOS  (X/SQRT  (PSQ  ) ) 

RETURN 

END 
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